{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "09dfe07e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sqlalchemy import create_engine, MetaData\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "from scipy import stats\n",
    "\n",
    "from shapely.geometry import Polygon\n",
    "from shapely.ops import transform\n",
    "import pyproj\n",
    "\n",
    "#from geoalchemy2 import Geometry  # <= not used but must be imported"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "62786daf",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3528: SAWarning: Skipped unsupported reflection of expression-based index unique_user_emails\n",
      "  util.warn(\n",
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3198: SAWarning: Did not recognize type 'geometry' of column 'extent'\n",
      "  util.warn(\n",
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3528: SAWarning: Skipped unsupported reflection of expression-based index unique_organization_names\n",
      "  util.warn(\n"
     ]
    }
   ],
   "source": [
    "api_engine = create_engine(f\"postgres://marxan-api:marxan-api@marxan-postgresql-api:5432/marxan-api\")\n",
    "api_meta = MetaData(schema=\"public\")\n",
    "api_meta.reflect(bind=api_engine, only=['output_results'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d9160ca7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3198: SAWarning: Did not recognize type 'geometry' of column 'the_geom'\n",
      "  util.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "dict_keys(['public.spatial_ref_sys', 'public.migrations', 'public.admin_regions_0', 'public.admin_regions', 'public.admin_regions_1', 'public.wdpa', 'public.admin_regions_2', 'public.features_data', 'public.planning_units_geom', 'public.planning_units_geom_square', 'public.planning_units_geom_hexagon', 'public.planning_units_geom_irregular', 'public.scenarios_pu_data', 'public.scenarios_pu_cost_data', 'public.output_results_data', 'public.scenario_features_data'])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geo_api_engine = create_engine(f\"postgres://marxan-geo-api:marxan-geo-api@marxan-postgresql-geo-api:5432/marxan-geo-api\")\n",
    "geo_api_meta = MetaData(schema=\"public\")\n",
    "geo_api_meta.reflect(bind=geo_api_engine)\n",
    "geo_api_meta.tables.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1937e48",
   "metadata": {},
   "source": [
    "[how many squares can we create in a certain extent?](https://www.redblobgames.com/grids)  \n",
    "\n",
    "#### Square grid maths:  \n",
    "65536 = 256*256 / 8 tile pixel extent per tiles in screen  \n",
    "grid = x km  \n",
    "extent = [min_x, min_y, max_x, max_y ]  \n",
    "width = max_x - min_x   \n",
    "height = max_y - min_y  \n",
    "n_rows = width / x   \n",
    "n_cols = height / x  \n",
    "ammount_cells = n_rows * n_cols  \n",
    "\n",
    "ammount_cells = extent_area / grid_area\n",
    "\n",
    "grid_area = Extent_area / ammount_cells  \n",
    "\n",
    "589824  = 65536 * 9 tiles in screen\n",
    "\n",
    "grid_area = Extent_area / 589824  \n",
    "**Difference seen in calculations related the real number are becouse in postgres/postgis we expect a regular indexed grid that only depends on gridsize. So as smaller our grid size better is our prediction. so we can use this functions to limit the ammount of geometries \n",
    "\n",
    "#### Hex grids maths:  \n",
    "\n",
    "w = sqrt(3) * size  \n",
    "h = 2 * size  \n",
    "The horizontal distance between adjacent hexagon centers is w.  \n",
    "The vertical distance between adjacent hexagon centers is h * 3/4\n",
    "\n",
    "\n",
    "\n",
    "#### Limitation of nº grid geometries to use in single computations; beyond those limints we should batch calculate stuff. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4798670f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "589824.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((256*256 / 1 ) * 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cae2045a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9216.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "589824.0/64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1bd7c4d1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9216.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((256*256 / 64 ) * 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bee10bde",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "170.46440972222223"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1571000/9216"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d957b0f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "b566a4be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5878547713829665"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zlevel = 7\n",
    "grid_size = 10\n",
    "m_per_pixel = 156412\n",
    "(math.sqrt(grid_size) * 1000) / (m_per_pixel/(2**zlevel))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "399f1acb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "78206/39103"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b54b4f09",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculateSquareOverExtent(geometry, grid_size_km):\n",
    "    \"\"\"Calculates the number of squares in grid covered by an extent and a grid_size \n",
    "    Parameters:\n",
    "    geometry (List(float)): [x_min, x_max, y_min, y_max] in epsg 4326\n",
    "    grid_size_km (float):  grid cell size in km\n",
    "\n",
    "    Returns:\n",
    "    int: Returning an aproach number of grids cells\n",
    "    \n",
    "    \"\"\"\n",
    "    grid_size_m = grid_size_km * 1000\n",
    "    wgs84 = pyproj.CRS('EPSG:4326')\n",
    "    webMercator = pyproj.CRS('EPSG:3857')\n",
    "    project = pyproj.Transformer.from_crs(wgs84, webMercator, always_xy=True).transform\n",
    "    utm_extent = transform(project, geometry).bounds\n",
    "    width = utm_extent[2] - utm_extent[0]\n",
    "    height = utm_extent[3] - utm_extent[1]\n",
    "\n",
    "    if width < grid_size_m or height < grid_size_m:\n",
    "        raise ValueError('grid_size is too big for extension')\n",
    "    nGrid = (round(width / grid_size_m) * round(height / grid_size_m))\n",
    "    areaExtent = width*height\n",
    "    areaGrid = grid_size_m ** 2\n",
    "    print(areaExtent)\n",
    "    return round(areaExtent/areaGrid)\n",
    "\n",
    "def calcMinGrid(extent_area):\n",
    "    return np.round(np.sqrt(extent_area / 589824))\n",
    "\n",
    "def calcMaxGrid(extent_area):\n",
    "    return np.round(np.sqrt(extent_area / 1))\n",
    "\n",
    "def calculateSquareOverExtent(geometry):\n",
    "    \"\"\"Calculates the number of squares in grid covered by an extent and a grid_size \n",
    "    Parameters:\n",
    "    geometry (List(float)): [x_min, x_max, y_min, y_max] in epsg 4326\n",
    "    grid_size_km (float):  grid cell size in km\n",
    "\n",
    "    Returns:\n",
    "    int: Returning an aproach number of grids cells\n",
    "    \n",
    "    \"\"\"\n",
    "    wgs84 = pyproj.CRS('EPSG:4326')\n",
    "    webMercator = pyproj.CRS('EPSG:3857')\n",
    "    project = pyproj.Transformer.from_crs(wgs84, webMercator, always_xy=True).transform\n",
    "    utm_extent = transform(project, geometry).bounds\n",
    "    width = utm_extent[2] - utm_extent[0]\n",
    "    height = utm_extent[3] - utm_extent[1]\n",
    "    areaExtent = width*height\n",
    "    min_grid = calcMinGrid(areaExtent)\n",
    "    max_grid = calcMaxGrid(areaExtent)\n",
    "    grid_size_m = round(min_grid/1000) * 1000\n",
    "    areaGrid = grid_size_m ** 2\n",
    "    print(f'geometry extent: {round(areaExtent/1000,2)} km2')\n",
    "    print(f'min grid size: {round(min_grid/1000)}')\n",
    "    print(f'max grid size: {round(max_grid/1000)}')\n",
    "    return round(areaExtent/areaGrid)\n",
    "\n",
    "def calculateHexOverExtent(geometry, grid_size_km):\n",
    "    \"\"\"Calculates the number of hexagons in grid covered by an extent and a grid_size \n",
    "    Parameters:\n",
    "    geometry (List(float)): [x_min, x_max, y_min, y_max] in epsg 4326\n",
    "    grid_size_km (float):  grid cell size in km\n",
    "\n",
    "    Returns:\n",
    "    int: Returning an aproach number of grids cells\n",
    "    \n",
    "    \"\"\"\n",
    "    grid_size_m = grid_size_km * 1000\n",
    "    wgs84 = pyproj.CRS('EPSG:4326')\n",
    "    webMercator = pyproj.CRS('EPSG:3857')\n",
    "    project = pyproj.Transformer.from_crs(wgs84, webMercator, always_xy=True).transform\n",
    "    utm_extent = transform(project, geometry).bounds\n",
    "    \n",
    "    width = utm_extent[2] - utm_extent[0]\n",
    "    height = utm_extent[3] - utm_extent[1]\n",
    "\n",
    "    if width < grid_size_m or height < grid_size_m:\n",
    "        raise ValueError('grid_size is too big for extension')\n",
    "    w = math.sqrt(3) * grid_size_m\n",
    "    h = 2 * grid_size_m\n",
    "    #The horizontal distance between adjacent hexagon centers is w. The vertical distance between adjacent hexagon centers is h * 3/4\n",
    "    return (round(width / w) * round(height / h))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ed42985d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def testQuery(polygon, size_km):\n",
    "    \"\"\"\n",
    "    \n",
    "    \"\"\"\n",
    "    size_m = size_km * 1000\n",
    "    geometry = f'{{\"type\":\"Polygon\",\"coordinates\":[{polygon}]}}'\n",
    "    query = f\"\"\"\n",
    "    select count(*) as count from \n",
    "    (SELECT (ST_SquareGrid({size_m}, ST_Transform(ST_GeomFromGeoJSON('{geometry}'), 3857))).* \n",
    "     ) grid\n",
    "    \"\"\"\n",
    "    with geo_api_engine.connect() as con:\n",
    "        rs = con.execute(query)\n",
    "        return rs.fetchall()[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e58bbd6a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "geometry extent: 380138877894.66 km2\n",
      "min grid size: 25\n",
      "max grid size: 19497\n",
      "608222\n",
      "609696\n",
      "2279\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"212.2453125\" height=\"127.26766756029625\" viewBox=\"-85.9078125 -44.1760626474805 212.2453125 127.26766756029625\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,38.91554226533524)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"2.0\" opacity=\"0.6\" d=\"M -78.046875,-36.3151251474805 L 118.47656249999999,-36.3151251474805 L 118.47656249999999,75.23066741281573 L -78.046875,75.23066741281573 L -78.046875,-36.3151251474805 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x7f98436d08e0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coord = [\n",
    "            [\n",
    "              -78.046875,\n",
    "              -36.3151251474805\n",
    "            ],\n",
    "            [\n",
    "              118.47656249999999,\n",
    "              -36.3151251474805\n",
    "            ],\n",
    "            [\n",
    "              118.47656249999999,\n",
    "              75.23066741281573\n",
    "            ],\n",
    "            [\n",
    "              -78.046875,\n",
    "              75.23066741281573\n",
    "            ],\n",
    "            [\n",
    "              -78.046875,\n",
    "              -36.3151251474805\n",
    "            ]\n",
    "          ]\n",
    "wgs84_pt = Polygon(coord)\n",
    "size_km = 25\n",
    "\n",
    "print(calculateSquareOverExtent(wgs84_pt))\n",
    "print(testQuery(coord, size_km))\n",
    "print(testQuery(coord, 416))\n",
    "wgs84_pt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "98b2ca7e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "640\n",
      "459\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-8.657666015625 41.277492002672126 4.4613281250000005 2.539717859639744\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,85.09470186498399)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"0.08922656250000001\" opacity=\"0.6\" d=\"M -7.03125,43.65197548731187 L -8.492431640625,41.44272637767212 L -4.361572265625,41.57436130598913 L -7.03125,43.65197548731187 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x7f98436bf8b0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wgs84_pt = Polygon([\n",
    "            [\n",
    "              -7.03125,\n",
    "              43.65197548731187\n",
    "            ],\n",
    "            [\n",
    "              -8.492431640625,\n",
    "              41.44272637767212\n",
    "            ],\n",
    "            [\n",
    "              -4.361572265625,\n",
    "              41.57436130598913\n",
    "            ],\n",
    "            [\n",
    "              -7.03125,\n",
    "              43.65197548731187\n",
    "            ]\n",
    "          ])\n",
    "\n",
    "query = \"\"\"\n",
    "select count(*) from \n",
    "(SELECT (ST_HexagonGrid(10000, ST_Transform(ST_GeomFromGeoJSON('{\"type\":\"Polygon\",\n",
    "        \"coordinates\":[[[-7.03125,43.65197548731187],[-8.492431640625,41.44272637767212],[-4.361572265625,41.57436130598913],[-7.03125,43.65197548731187]]]}'), 3857))).* \n",
    " ) grid\n",
    "\"\"\"\n",
    "with geo_api_engine.connect() as con:\n",
    "    rs = con.execute(query)\n",
    "    for row in rs:\n",
    "        print(row[0])\n",
    "\n",
    "print(calculateHexOverExtent(wgs84_pt, 10))\n",
    "wgs84_pt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "44bf8672",
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "calculateSquareOverExtent() takes 1 positional argument but 2 were given",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-8-2a44ff0a9f2d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# 100 linearly spaced numbers represent km\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mx0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0my0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcalculateSquareOverExtent\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwgs84_pt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0my1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mtestQuery\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoord\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-8-2a44ff0a9f2d>\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# 100 linearly spaced numbers represent km\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mx0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0my0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcalculateSquareOverExtent\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwgs84_pt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0my1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mtestQuery\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoord\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: calculateSquareOverExtent() takes 1 positional argument but 2 were given"
     ]
    }
   ],
   "source": [
    "coord = [[-7.03125,43.65197548731187],[-8.492431640625,41.44272637767212],[-4.361572265625,41.57436130598913],[-7.03125,43.65197548731187]]\n",
    "wgs84_pt = Polygon(coord)\n",
    "# 100 linearly spaced numbers represent km\n",
    "x0 = np.linspace(100,1,100)\n",
    "y0 = [calculateSquareOverExtent(wgs84_pt, x) for x in x0]\n",
    "y1 = [testQuery(coord, x) for x in x0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a5bbddfa",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'y0' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-24174dbc0997>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;31m# plot the function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"line 1\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     14\u001b[0m \u001b[0;31m# plot the function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'g'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"line 2\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'y0' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAATOUlEQVR4nO3df4xlZ13H8feH3Zbyoz/panBbaMGCrqaFdqgoCC2otFVTQSQtlcb6Y7NKpQlRC/4mxIjBktpQrCuuFBKpRhtYTKXBH1C1VnfQpe1uLVkLdMeSdEuR0hap2379495lLtPZfe6ZnXPn7uz7ldxkzj3PPfd7n8zcz5znnPOcVBWSJB3IU1a6AEnS9DMsJElNhoUkqcmwkCQ1GRaSpCbDQpLU1GtYJNmS5P4kd+5nfZJck2RXktuTnNlnPZKkpel7z+IDwHkHWH8+cNrwsRH4w57rkSQtQa9hUVW3AA8eoMmFwAdr4DbguCTP7rMmSVJ3a1f4/dcDu0eW54bPfXFhwyQbGex9sGHDhrN27NgxkQIlaRXJUl+40ge4Fyt80flHqmpzVc1U1czTnva0nsuSJI1a6bCYA04eWT4JuG+FapEk7cdKh8VW4NLhWVEvBb5SVU8agpIkraxej1kk+TBwDnBikjngt4AjAKrqOuAm4AJgF/AocFmf9UiSlqbXsKiqixvrC3hznzVIkg7eSg9DSZIOAYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDX1HhZJzktyd5JdSd62yPpjk3wsyWeS7EhyWd81SZK66TUskqwBrgXOBzYAFyfZsKDZm4GdVXUGcA5wVZIj+6xLktRN33sWZwO7quqeqnoMuAG4cEGbAo5OEuCZwIPA3p7rkiR10HdYrAd2jyzPDZ8b9V7gO4H7gDuAK6rqiYUbSrIxyWyS2T179vRVryRpEX2HRRZ5rhYsvwbYDnwb8CLgvUmOedKLqjZX1UxVzaxbt26565QkHUDfYTEHnDyyfBKDPYhRlwE31sAu4HPAd/RclySpg77DYhtwWpJThwetLwK2LmhzL/BqgCTfCrwQuKfnuiRJHaztc+NVtTfJ5cDNwBpgS1XtSLJpuP464J3AB5LcwWDY6sqqeqDPuiRJ3aRq4SGE6TczM1Ozs7MrXYYkHWoWO448Fq/gliQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqWlJYZHkKUmOWe5iJEnTaeywSPJnSY5J8gxgJ3B3kl/urzRJ0rTosmexoaoeAn4MuAl4DvCmPoqSJE2XLmFxRJIjGITFR6vq/4DqpSpJ0lTpEhZ/BHweeAZwS5LnAg/1UZQkabqsHbdhVV0DXDPy1BeSnLv8JUmSpk0zLJK8tdHkPctUiyRpSo2zZ3F071VIkqZaMyyq6h2TKESSNL3GGYa65kDrq+oty1eOJGkajTMM9emDeYMk5wF/AKwB3l9V71qkzTnA1cARwANV9cqDeU9J0vIaZxjq+tHlJM+oqkfG2XiSNcC1wA8Cc8C2JFuraudIm+OA9wHnVdW9Sb6lQ/2SpAnoMt3H9ybZCdw1XD4jyfsaLzsb2FVV91TVY8ANwIUL2rwRuLGq7gWoqvvHrl6SNBFdLsq7GngN8CWAqvoM8IrGa9YDu0eW54bPjXoBcHySTyb5dJJLF9tQko1JZpPM7tmzp0PZkqSD1WnW2araveCpxxsvyWKbWbC8FjgL+GEGYfQbSV6wyHtvrqqZqppZt27duCVLkpbB2FdwA7uTfB9QSY4E3sJwSOoA5oCTR5ZPAu5bpM0Dw+MgjyS5BTgD+GyH2iRJPeqyZ7EJeDODYaQ54EXD5QPZBpyW5NRhwFwEbF3Q5qPA9ydZm+TpwPfQDiFJ0gR1mRvqAeCSLhuvqr1JLgduZnDq7Jaq2pFk03D9dVV1V5KPA7cDTzA4vfbOLu8jSepXqsabZTzJ9cAVVfU/w+Xjgauq6qf7K29xMzMzNTs7O+m3laRD3WLHkcfSZRjq9H1BAVBVXwZevNQ3liQdOrqExVOGexMAJDmBbgfIJUmHqC5f9lcBtyb5Swanv74B+J1eqpIkTZUuB7g/mGQWeBWDca/XLZi24/jh0JQkaZXpNIw0DIed+1n9d8CZB12RJGnqdLqCu2HJR9klSdNtOcNivHNwJUmHnOUMC0nSKuUwlCSpqcv9LD7UeO7Vy1KRJGnqdNmz+K7RheFd8M7at1xVDy5XUZKk6dIMiyRvT/JV4PQkDw0fXwXuZzBjrCRplWuGRVX9blUdDby7qo4ZPo6uqmdV1dsnUKMkaYV1uYL77UnWA88dfV1V3dJHYZKk6TF2WCR5F4ObF+1k/naqBRgWkrTKdZnu47XAC6vq630VI0maTl3OhroHOKKvQiRJ06vLnsWjwPYkfwd8Y++iqt6y7FVJkqZKl7DYOnxIkg4zXc6Gur7PQiRJ06vL2VCfY5GZZavqectakSRp6nQZhpoZ+fko4CeAE5a3HEnSNBr7bKiq+tLI47+r6moGt1iVJK1yXYahRm+Z+hQGexpHL3tFkqSp02UY6qqRn/cCnwfesKzVSJKmUpezoc7tsxBJ0vTqcvOjY5O8J8ns8HFVkmP7LE6SNB26TPexBfgqg6GnNwAPAX/aR1GSpOnS5ZjF86vqx0eW35Fk+zLXI0maQl32LL6W5OX7FpK8DPja8pckSZo2XfYsfh64fnicIsCDwE/1UZQkabp0ORtqO3BGkmOGyw/1VZQkabp0uSjvOOBS4BRgbRLAKcol6XDQZRjqJuA24A7giX7KkSRNoy5hcVRVvbW3SiRJU6vL2VAfSvJzSZ6d5IR9j9aLkpyX5O4ku5K87QDtXpLk8SSv71CTJGkCuuxZPAa8G/g15u9rUcB+72eRZA1wLfCDwBywLcnWqtq5SLvfA27uUI8kaUK6hMVbgW+vqgc6vOZsYFdV3QOQ5AbgQmDngna/CPwV8JIO25YkTUiXYagdwKMdt78e2D2yPDd87huSrAdeC1x3oA0l2bhvXqo9e/Z0LEOSdDC67Fk8DmxP8g/A1/c92Th1Nos8t/DWrFcDV1bV4/tOx11MVW0GNgPMzMw86faukqT+dAmLjwwfXcwBJ48snwTct6DNDHDDMChOBC5Isrequr6XJKknXa7gvn4J298GnJbkVOC/gYuANy7Y7qn7fk7yAeCvDQpJmi5d7mfxI0n+I8mDSR5K8tUkB5zyo6r2ApczOMvpLuAvqmpHkk1JNh1c6ZKkSUnVeMP/SXYBrwPuqHFf1JOZmZmanZ1dyRIk6VC0/wPDDV3OhtoN3LnSQSFJmrwuB7h/Bbgpyaf45rOh3rPsVUmSpkqXsPgd4GHgKODIfsqRJE2jLmFxQlX9UG+VSJKmVpdjFn+bxLCQpMNQl7B4M/DxJF8b99RZSdLq0OWivKOHU5KfxuC4hSTpMNHltqo/C1zBYMqO7cBLgVuBV/dSmSRpanQZhrqCwRTiX6iqc4EXA12mK5ckHaK6hMX/VtX/AiR5alX9J/DCfsqSJE2TLqfOziU5jsHMs59I8mWePIOsJGkV6nKA+7XDH397eE+LY4GP91KVJGmqdNmz+Iaq+tRyFyJJml5djllIkg5ThoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktTUe1gkOS/J3Ul2JXnbIusvSXL78HFrkjP6rkmS1E2vYZFkDXAtcD6wAbg4yYYFzT4HvLKqTgfeCWzusyZJUnd971mcDeyqqnuq6jHgBuDC0QZVdWtVfXm4eBtwUs81SZI66jss1gO7R5bnhs/tz88Af7PYiiQbk8wmmd2zZ88ylihJauk7LLLIc7Vow+RcBmFx5WLrq2pzVc1U1cy6deuWsURJUsvanrc/B5w8snwScN/CRklOB94PnF9VX+q5JklSR33vWWwDTktyapIjgYuAraMNkjwHuBF4U1V9tud6JElL0OueRVXtTXI5cDOwBthSVTuSbBquvw74TeBZwPuSAOytqpk+65IkdZOqRQ8hTLWZmZmanZ1d6TIk6VCz2HHksXgFtySpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpKbewyLJeUnuTrIrydsWWZ8k1wzX357kzL5rkiR102tYJFkDXAucD2wALk6yYUGz84HTho+NwB/2WZMkqbu+9yzOBnZV1T1V9RhwA3DhgjYXAh+sgduA45I8u+e6JEkdrO15++uB3SPLc8D3jNFmPfDF0UZJNjLY8wD4epI7l7fUQ9aJwAMrXcSUsC/m2Rfz7It5d1bVdy/lhX2HRRZ5rpbQhqraDGwGSDJbVTMHX96hz76YZ1/Msy/m2Rfzkswu9bV9D0PNASePLJ8E3LeENpKkFdR3WGwDTktyapIjgYuArQvabAUuHZ4V9VLgK1X1xYUbkiStnF6Hoapqb5LLgZuBNcCWqtqRZNNw/XXATcAFwC7gUeCyMTa9uaeSD0X2xTz7Yp59Mc++mLfkvkjVkw4PSJL0TbyCW5LUZFhIkpqmOiycKmTeGH1xybAPbk9ya5IzVqLOSWj1xUi7lyR5PMnrJ1nfJI3TF0nOSbI9yY4kn5p0jZMyxt/IsUk+luQzw74Y5/joISfJliT37+9atCV/b1bVVD4YHBD/L+B5wJHAZ4ANC9pcAPwNg2s1Xgr860rXvYJ98X3A8cOfzz+c+2Kk3d8zOIHi9Std9wr+XhwH7ASeM1z+lpWuewX74leB3xv+vA54EDhypWvvoS9eAZzJ4AK8xdYv6XtzmvcsnCpkXrMvqurWqvrycPE2BterrEbj/F4A/CLwV8D9kyxuwsbpizcCN1bVvQBVtVr7Y5y+KODoJAGeySAs9k62zP5V1S0MPtv+LOl7c5rDYn/TgHRtsxp0/Zw/w+A/h9Wo2RdJ1gOvBa6bYF0rYZzfixcAxyf5ZJJPJ7l0YtVN1jh98V7gOxlc9HsHcEVVPTGZ8qbKkr43+57u42As21Qhq8DYnzPJuQzC4uW9VrRyxumLq4Erq+rxwT+Rq9Y4fbEWOAt4NfA04F+S3FZVn+27uAkbpy9eA2wHXgU8H/hEkn+sqod6rm3aLOl7c5rDwqlC5o31OZOcDrwfOL+qvjSh2iZtnL6YAW4YBsWJwAVJ9lbVRyZS4eSM+zfyQFU9AjyS5BbgDGC1hcU4fXEZ8K4aDNzvSvI54DuAf5tMiVNjSd+b0zwM5VQh85p9keQ5wI3Am1bhf42jmn1RVadW1SlVdQrwl8AvrMKggPH+Rj4KfH+StUmezmDW57smXOckjNMX9zLYwyLJtwIvBO6ZaJXTYUnfm1O7Z1H9TRVyyBmzL34TeBbwvuF/1HtrFc60OWZfHBbG6YuquivJx4HbgSeA91fVqpvef8zfi3cCH0hyB4OhmCuratVNXZ7kw8A5wIlJ5oDfAo6Ag/vedLoPSVLTNA9DSZKmhGEhSWoyLCRJTYaFJKnJsJAkNRkW0hiSbFpsqowkp+xvds/9bOemJMcta3HSBHjqrNSQZG1VLTrhXJJTgL+uqu+ebFXSZE3tRXnSJCT5DeASBhOrPQB8uqp+P8kngVuBlwFbkxwNPDxcdxawhcEFTf+0n+0+G/hz4BgGf2c/X1X/mOTzDKYjeT2wadj8WODzVXVukh8C3gE8lcGU25dV1cPL/8mlbhyG0mEryQzw48CLgdcx+BIfdVxVvbKqrlrw/J8Cb6mq7z3A5t8I3FxVL2IwF9P20ZXDq6tfBLyEwVw970lyIvDrwA9U1ZnALPDWJXw0adm5Z6HD2cuBj1bV1wCSfGzB+j9f+IIkxzIIkX13nPsQg5tNLbQN2JLkCOAjVbV9PzX8AfD3VfWxJD8CbAD+eThly5HAv3T7SFI/DAsdzlrzlz+yn9c0D/RV1S1JXgH8MPChJO+uqg9+04aSnwKeC1w+su1PVNXFre1Lk+YwlA5n/wT8aJKjkjyTwRf7AVXV/wBfSbLvfiGXLNYuyXOB+6vqj4E/YXCby9H1ZwG/BPzkyA14bgNeluTbh22enuQF3T+WtPzcs9Bhq6q2JdnK4H7NX2BwjOArY7z0MgZDTI8ymOV0MecAv5zk/4CHgYWn3V4OnAD8w3DIabaqfna4t/HhJE8dtvt1Vt+9J3QI8tRZHdaSPLOqHh7e6+EWYGNV/ftK1yVNG/csdLjbnGQDcBRwvUEhLc49C0lSkwe4JUlNhoUkqcmwkCQ1GRaSpCbDQpLU9P/1x4wS1oiUiwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.xlabel(\"grid size\")\n",
    "plt.ylabel(\"ammount_cells\")\n",
    "\n",
    "# plot the function\n",
    "plt.plot(x0, y0, 'r', label = \"line 1\")\n",
    "# plot the function\n",
    "plt.plot(x0, y1, 'g', label = \"line 2\")\n",
    "\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1cf41152",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'y1' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-9807fbeeee47>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;31m# plot the function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;31m# show the plot\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'y1' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAATrklEQVR4nO3df5BlZX3n8ffHGQkQFVgYjM4MkZBRM5sChRZRE0WziQyaYjWagJaUSGrCrqjJblJQqYquMalN3HLLFUEyIQimsmIirCGUwqZ2y5ASydK68mMgmMkYYQIVhkgwCylh4Lt/nAt9bXqePtMz5/ad6ferqqvvuee55377qe776XOec56TqkKSpN151nIXIEmabgaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaBguKJJcneSDJHbtZnySfSLItyW1JThyqFknS0g25R3EFcFpj/SZgw+hrM/CpAWuRJC3RYEFRVTcC32k0OQP4THVuBg5P8oKh6pEkLc3qZXzvtcC9Y8s7Rs/dP79hks10ex1s3LjxpK1bt06kQEk6gGSpL1zOweyFil5wPpGq2lJVM1U1c8ghhwxcliRp3HIGxQ5g/djyOuC+ZapFkrQbyxkU1wJnj85+OgV4uKqecdhJkrS8BhujSPJZ4FTgqCQ7gA8BzwaoqkuBLwKnA9uAR4FzhqpFkrR0gwVFVZ21yPoC3jvU+0uS9g2vzJYkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktQ0aFAkOS3J3Um2JblwgfWHJfmzJLcm2ZrknCHrkSTtucGCIskq4GJgE7AROCvJxnnN3gvcWVUnAKcCH0ty0FA1SZL23JB7FCcD26pqe1U9BlwFnDGvTQHPTRLgOcB3gF0D1iRJ2kNDBsVa4N6x5R2j58Z9Evgx4D7gduADVfXk/A0l2ZxkNsnszp07h6pXkrSAIYMiCzxX85bfCHwDeCHwMuCTSZ73jBdVbamqmaqaWbNmzb6uU5LUMGRQ7ADWjy2vo9tzGHcOcE11tgHfAl46YE2SpD00ZFDcAmxIcuxogPpM4Np5be4BfgogyfOBlwDbB6xJkrSHVg+14araleR84AZgFXB5VW1Nct5o/aXAR4ArktxOd6jqgqp6cKiaJEl7LlXzhw2m28zMTM3Ozi53GZK0v1lo3LgXr8yWJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaFg2KJIcm+Y0kvz9a3pDkzcOXJkmaBn32KD4NfA941Wh5B/Bbg1UkSZoqfYLiuKr6KPA4QFX9C5BBq5IkTY0+QfFYkkOAAkhyHN0ehiRpBVjdo82HgOuB9Un+CHgN8O4hi5IkTY9Fg6Kq/jzJ14FT6A45faCqHhy8MknSVFg0KJK8dvTwn0ffNyahqm4crixJ0rToc+jp18YeHwycDHwNeMMgFUmSpkqfQ08/O76cZD3w0cEqkiRNlaVcmb0D+PF9XYgkaTr1GaO4iNGpsXTB8jLg1gFrkiRNkT5jFLNjj3cBn62qrwxUjyRpyvQZo7hyEoVIkqbTboMiye3MHXL6vlVAVdXxg1UlSZoarT2KvZ4hNslpwH8DVgGXVdXvLNDmVODjwLOBB6vqdXv7vpKkfWe3QVFV396bDSdZBVwM/DTdmVK3JLm2qu4ca3M4cAlwWlXdk+TovXlPSdK+1+d+FKckuSXJ/0vyWJInkny3x7ZPBrZV1faqegy4CjhjXpt3ANdU1T0AVfXAnv4AkqRh9bmO4pPAWcDfAIcAvwhc1ON1a4F7x5Z3jJ4b92LgiCRfTvK1JGcvtKEkm5PMJpnduXNnj7eWJO0rvS64q6ptwKqqeqKqPg28vsfLFrpnxfzB8dXAScCbgDcCv5HkxQu8/5aqmqmqmTVr1vQpWZK0j/S5juLRJAcB30jyUeB+4Ad7vG4HsH5seR1w3wJtHqyqR4BHktwInAB8s8f2JUkT0GeP4l2jducDj9B9+P9cj9fdAmxIcuwoaM4Erp3X5k+Bn0yyOsmhwCuBu/oWL0kaXp89ihOBL1bVd4EP991wVe1Kcj5wA93psZdX1dYk543WX1pVdyW5HrgNeJLuFNo79vinkCQNJlULXVM31iD5NN2U4jfSnbl0Q1XtmkBtC5qZmanZ2dnFG0qSxi00btzLooeequoc4EeBP6E7nfVvk1y21DeUJO1f+hx6oqoeT/IlurOWDqG7HuIXhyxMkjQd+lxwd1qSK4BtwNuAy4AXDFyXJGlK9NmjeDfd2MQvVdX3hi1HkjRt+kwzfuYkCpEkTael3ApVkrSCGBSSpCaDQpLUtJQ73AHgHe4kaWXoc4e7946+/+Ho+zuBRwerSJI0VRa9w12S11TVa8ZWXZjkK8BvDl2cJGn59Rmj+MEkP/HUQpJX02+acUnSAaDPBXfnApcnOYxuzOJh4D2DViVJmhp9Lrj7GnBCkufRzTb78PBlSZKmRZ+5np6f5A+Az1XVw0k2Jjl3ArVJkqZAnzGKK+huPvTC0fI3gV8eqB5J0pTpExRHVdUf092BjtFNi54YtCpJ0tToExSPJDmS0cV3SU6hG9CWJK0Afc56+g/AtcBxo+sn1gBvH7QqSdLU6BMUW4HXAS+hu+fq3ThHlCStGH0+8L9aVbuqamtV3VFVjwNfHbowSdJ0aE0K+EPAWuCQJC+n25sAeB5w6ARqkyRNgdahpzfS3QZ1HfBfx57/Z+DXB6xJkjRFWpMCXglcmeTnqurqCdYkSZoifabwuDrJm4B/DRw89ryzx0rSCtBnCo9LgV8A3kc3TvF24IcHrkuSNCX6nPX06qo6G3ioqj4MvApYP2xZkqRp0Sco/mX0/dEkLwQeB44driRJ0jTpc8HddUkOB/4L8HW6qTwuG7IoSdL06DOY/ZHRw6uTXAcc7D0pJGnlaF1w99bGOqrqmmFKkiRNk9Yexc821hVgUEjSCtC64O6cSRYiSZpOi45RJPngQs97wZ0krQx9znp6ZOzxwcCbgbuGKUeSNG0WvY6iqj429vXbwKl0s8ouKslpSe5Osi3JhY12r0jyRJK39a5ckjQRS7kB0aHAjyzWKMkq4GJgE7AROCvJxt20+13ghiXUIkkaWJ8xitsZ3S8bWEV3K9Q+4xMnA9uqavtoO1cBZwB3zmv3PuBq4BU9a5YkTVCfMYo3jz3eBfxDVe3q8bq1wL1jyzuAV443SLIWeAvwBhpBkWQzsBngmGOO6fHWkqR9pc+hpxcA36mqb1fV3wMHJ3nlYi9i7o5442re8seBC6rqidaGqmpLVc1U1cyaNWt6vLUkaV/ps0fxKeDEseVHF3huITv4/llm1wH3zWszA1yVBOAo4PQku6rqCz3qkiRNQJ+gSFU9vSdQVU8m6fO6W4ANSY4F/h44E3jHeIOqenoW2iRXANcZEpI0Xfocetqe5P1Jnj36+gCwfbEXjcYxzqc7m+ku4I+ramuS85Kct3dlS5ImJWM7Cws3SI4GPkE34FzA/wJ+uaoeGL68Z5qZmanZ2dnleGtJ2p8tNG7cS59pxh+gO2wkSVqBWtOMX8Qzz1J6WlW9f5CKJElTpbVH4fEdSVJzmvErJ1mIJGk69ZnCYw1wAd18TQc/9XxVvWHAuiRJU6LP6bF/RHd667HAh4G/o7tGQpK0AvQJiiOr6g+Ax6vqL6rqPcApA9clSZoSfa6wfnz0/f4kb6KbhmPdcCVJkqZJn6D4rSSHAf8RuAh4HvArg1YlSZoafS64u2708GHg9cOWI0maNouOUSS5MsnhY8tHJLl80KokSVOjz2D28VX1T08tVNVDwMsHq0iSNFX6BMWzkhzx1EKSf0W/sQ1J0gGgzwf+x4Cbknyebu6nnwd+e9CqJElTo89g9meSzNJNMx7grVV15+CVSZKmQq9DSKNgMBwkaQXqM0YhSVrBDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklS06BBkeS0JHcn2ZbkwgXWvzPJbaOvm5KcMGQ9kqQ9N1hQJFkFXAxsAjYCZyXZOK/Zt4DXVdXxwEeALUPVI0lamiH3KE4GtlXV9qp6DLgKOGO8QVXdVFUPjRZvBtYNWI8kaQmGDIq1wL1jyztGz+3OucCXFlqRZHOS2SSzO3fu3IclSpIWM2RQZIHnasGGyevpguKChdZX1ZaqmqmqmTVr1uzDEiVJi1k94LZ3AOvHltcB981vlOR44DJgU1X944D1SJKWYMg9iluADUmOTXIQcCZw7XiDJMcA1wDvqqpvDliLJGmJBtujqKpdSc4HbgBWAZdX1dYk543WXwp8EDgSuCQJwK6qmhmqJknSnkvVgsMGU2tmZqZmZ2eXuwxJ2t8sNG7ci1dmS5KaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJaho0KJKcluTuJNuSXLjA+iT5xGj9bUlOHLIeSdKeGywokqwCLgY2ARuBs5JsnNdsE7Bh9LUZ+NRQ9UiSlmbIPYqTgW1Vtb2qHgOuAs6Y1+YM4DPVuRk4PMkLBqxJkrSHVg+47bXAvWPLO4BX9mizFrh/vFGSzXR7HADfS3LHvi11v3UU8OByFzEl7Is59sUc+2LOHVX140t54ZBBkQWeqyW0oaq2AFsAksxW1czel7f/sy/m2Bdz7Is59sWcJLNLfe2Qh552AOvHltcB9y2hjSRpGQ0ZFLcAG5Icm+Qg4Ezg2nltrgXOHp39dArwcFXdP39DkqTlM9ihp6raleR84AZgFXB5VW1Nct5o/aXAF4HTgW3Ao8A5PTa9ZaCS90f2xRz7Yo59Mce+mLPkvkjVM4YEJEl6mldmS5KaDApJUtPUBoXTf8zp0RfvHPXBbUluSnLCctQ5CYv1xVi7VyR5IsnbJlnfJPXpiySnJvlGkq1J/mLSNU5Kj7+Rw5L8WZJbR33RZzx0v5Pk8iQP7O5asyV/blbV1H3RDX7/LfAjwEHArcDGeW1OB75Edy3GKcBfLXfdy9gXrwaOGD3etJL7Yqzd/6Y7WeJty133Mv5eHA7cCRwzWj56uetexr74deB3R4/XAN8BDlru2gfoi9cCJ9JdXLfQ+iV9bk7rHoXTf8xZtC+q6qaqemi0eDPd9SgHoj6/FwDvA64GHphkcRPWpy/eAVxTVfcAVNWB2h99+qKA5yYJ8By6oNg12TKHV1U30v1su7Okz81pDYrdTe2xp20OBHv6c55L9x/DgWjRvkiyFngLcOkE61oOfX4vXgwckeTLSb6W5OyJVTdZffrik8CP0V3Qezvwgap6cjLlTZUlfW4OOYXH3thn038cAHr/nEleTxcUPzFoRcunT198HLigqp7o/nk8YPXpi9XAScBPAYcAX01yc1V9c+jiJqxPX7wR+AbwBuA44M+T/GVVfXfg2qbNkj43pzUonP5jTq+fM8nxwGXApqr6xwnVNml9+mIGuGoUEkcBpyfZVVVfmEiFk9P3b+TBqnoEeCTJjcAJwIEWFH364hzgd6o7UL8tybeAlwL/ZzIlTo0lfW5O66Enp/+Ys2hfJDkGuAZ41wH43+K4Rfuiqo6tqhdV1YuAzwP//gAMCej3N/KnwE8mWZ3kULrZm++acJ2T0Kcv7qHbsyLJ84GXANsnWuV0WNLn5lTuUdRw03/sd3r2xQeBI4FLRv9J76oDcMbMnn2xIvTpi6q6K8n1wG3Ak8BlVXXATdHf8/fiI8AVSW6nO/xyQVUdcNOPJ/kscCpwVJIdwIeAZ8PefW46hYckqWlaDz1JkqaEQSFJajIoJElNBoUkqcmgkCQ1GRTSHkjyd0mOmpbtSJNgUGhFGl1w5O+/1IN/KFoxkrwoyV1JLgG+DqxP8mtJbhnNzf/hsbZfGE2ktzXJ5kW2+++SfHRs+d1JLuqznVFNd4wt/2qS/zR6fFyS60ev/8skL93rTpCWwKDQSvMSummWXz56vIFumuqXASclee2o3Xuq6iS6uaPen+TIxjY/D7x1bPkXgM8tYTvzbQHeN3r9rwKX7MFrpX1mKqfwkAb07dE8/AA/M/r6v6Pl59AFx410H+pvGT2/fvT8gpMtVtXOJNtHc+f8DV0AfWW0uvd2xiV5Dt0Nqf5kbBbcH+j1E0r7mEGhleaRsccB/nNV/d54gySnAv8GeFVVPZrky8DBi2z3c8DPA38N/I+qqp7b2cX379k/tf5ZwD9V1cv6/FDSkDz0pJXsBuA9o//eSbI2ydHAYcBDow/3l9LdMnIx1wD/FjiLucNOfbbzD8DRSY5M8gPAmwFG90n4VpK3j2pLDuB7oWu6GRRasarqfwL/ne6GPrfTjTU8F7geWJ3kNrpZR2/e/Vae3tZDdPen/uGqeuoeB4tup6oeB34T+CvgOro9kqe8Ezg3ya3AVha+7as0OGePlSQ1uUchSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKa/j+pWpL/tXqGlwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.xlabel(\"real value\")\n",
    "plt.ylabel(\"calculated value\")\n",
    "\n",
    "# plot the function\n",
    "plt.plot(y1, y0, 'r')\n",
    "\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5c3c7916",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'y0' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-21f8c9068251>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0merror\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mfig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_subplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mspines\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'left'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_position\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'zero'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'y0' is not defined"
     ]
    }
   ],
   "source": [
    "error = ((np.asarray(y0) - np.asarray(y1))/np.asarray(y1))*100\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.xlabel(\" Extent / km2 grid size km \")\n",
    "plt.ylabel(\"% error between calculated and real\")\n",
    "\n",
    "# plot the function\n",
    "plt.plot(153531695409.90253/x0, error, 'r')\n",
    "\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "c4491981",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApFklEQVR4nO3deZgU5bn38e8NDIsICIKEsRWRQQRkUQcU5DXBRAGNoolGkhgggiRED8aYRbO8ZjnJMceTvOrJkTjGKMSFgxsQxVEkRo2KOKO4AJLBgDIDyqaggCjD/f7x1EAzztKzdNfM9O9zXXVV91NV3Xc3Rd9T9Wzm7oiISPZqFXcAIiISLyUCEZEsp0QgIpLllAhERLKcEoGISJZrE3cA9VTvpk7jxo2jsLCwMWMREWkurKrCrLsi2LJlS9whiIg0KVmXCERE5GBKBCIiWU6JQEQkyykRiIhkOSUCEZEsl7ZEYGZ/NrNNZvZ6Ulk3M1tsZiXRumvStmvNbI2ZrTazsUnlJ5vZa9G2m83MAPbs2cPFF19MXl4ep5xyCuvWrUvXRxERadHSeUVwJzCuUtk1wBJ37wcsiZ5jZgOBicCg6JhbzKx1dMwsYDrQL1rGAdx+++107dqVNWvWcNVVV/GjH/0ojR9FRKTlSlsicPengW2ViicAs6PHs4Hzk8rnuvsed18LrAFGmFkvoLO7P+9hvOw5FccsWLCAyZMnA3DhhReyZMkSah1S++abobS0gZ9MRKRlyXQdQU933wgQrY+Iyo8E1iftVxqVHRk9rlxOWVkZRx11FABt2rShS5cubN26tco3LSgoID8/n8d//nPKq9lHRCRbNZXK4qq6PXsN5VX+9R9VH3zK9OnTKSoq4qxzzqF1q6bykUVEmoZM/yq+G93uIVpvispLgaOS9ksAG6LyRBXlJBIJ1q8PFxF79+5l+/btdOvWreZ3z8kBzcgmInKQTCeChcDk6PFkYEFS+UQza2dmfQiVwsui20cfmNmpUWuhSRXHnHfeecyeHaob7r//fs4444xqrwj2UyIQEfmUtI0+amb3Ap8DuptZKXAdcD0wz8ymAm8DFwG4+wozmwesBPYCl7t7efRSMwgtkDoAj0YLU6dO5Rvf+AZ5eXl069aNuXPn1h5Uu3ZKBCIilVgznby+fkH/6Efk33ADRfv2NXI4IiLNgoahpkOHcEWgRCAisl/2JQKAjz6KNw4RkSYkuxJBx45h/eGH8cYhItKEZFci6NQprD/4IN44RESaECUCEZEsl12JoHPnsFYiEBHZL7sSwWGHhfX778cZhYhIk5JdiaBrNP3BtsqDooqIZK/sSgQVYxG99168cYiINCHZlQi6dAlrDUUtIrJfdiWCVq2gTRvYsiXuSEREmozsSgQQEsGmTbXvJyKSJZQIRESyXPYlgpwcePfduKMQEWkysjMRbNwYdxQiIk1GdiaCDz+EHTvijkREpEnIvkTQtm1Yb9gQbxwiIk1E9iWCnJywLi2NNw4RkSYi+xJBu3Zh/dZb8cYhItJEZF8iyMkJHcvefjvuSEREmoTsSwRmkJurKwIRkUj2JQKAPn3gX/+KOwoRkSYhOxNB376wZk3cUYiINAnZmQjy8kKnsp07445ERCR22ZsIAN58M944RESagOxMBP37h/Xq1fHGISLSBGRnIjjuuNB66I034o5ERCR22ZkIDjkEeveGVavijkREJHbZmQgABgyAlSvjjkJEJHbZmwgGDw5XBJ98EnckIiKxyu5E8PHHUFISdyQiIrGKJRGY2VVmtsLMXjeze82svZl1M7PFZlYSrbsm7X+tma0xs9VmNraivLi4mMGDB5OXl8fMmTNx99SDOOGEsH711cb7YCIizVDGE4GZHQnMBPLd/QSgNTARuAZY4u79gCXRc8xsYLR9EDAOuKW8vByAGTNmUFBQQElJCSUlJRQWFqYeyMCBYQC65csb66OJiDRLcd0aagN0MLM2wCHABmACMDvaPhs4P3o8AZjr7nvcfS2wZtmyZWzcuJEdO3YwcuRIzIxJkyYxf/781CNo2zbcHnrppcb5RCIizVTGE4G7lwH/BbwNbAS2u/vjQE933xjtsxE4IjrkSGB90kuUlpWVUVZWRiKR2F+YSCQoKyur8j0LCgrIz88nPz+fzZs3H9hw0kkhEdTllpKISAsTx62hroS/8vsAuUBHM7ukpkOqeI0q6wPMPrUrANOnT6eoqIiioiJ69OhxYMNJJ8HWrRqSWkSyWhy3hr4ArHX3ze7+CfAgMAp418x6AUTrTdH+pcBRSccncnNzSSQSlCZNN1laWkpubm7dIjnllLB+4YV6fRARkZYgjkTwNnCqmR1i4U/4zwOrgIXA5GifycCC6PFCYKKZtTOzPkC/ESNG0KtXLzp16sTSpUtxd+bMmcOECRPqFsngwdC+vRKBiGS1Npl+Q3d/wczuB14C9gIvAwXAocA8M5tKSBYXRfuvMLN5wMpo/8tbt269CGDWrFlMmTKF3bt3M378eMaPH1+3YHJywu2hpUsb6dOJiDQ/Vqe2901HvYPOz8+nqKjoQMEPfgA33wzbt4erAxGRlqvKitTs7VlcYfTo0MM4OTmIiGQRJYLTTgvrZ56JNw4RkZgoEXTvHnoZP/103JGIiMRCiQBgzJhwRaCRSEUkCykRQEgEO3fCiy/GHYmISMYpEQB89rNh/be/xRuHiEgMlAgg1BOcdBI8/njckYiIZJwSQYWxY+H552HHjrgjERHJKCWCCmedBXv36vaQiGQdJYIKo0ZBp06waFHckYiIZJQSQYW2bcPtoYcf1vwEIpJVlAiSffGLsHEjvPxy3JGIiGSMEkGys8+GVq1gwYLa9xURaSGUCJL16BEGoXvoobgjERHJGCWCyi64AF57DUpK4o5ERCQjlAgqu+CCsH7ggXjjEBHJECWCynr3DnMZz5sXdyQiIhmhRFCViRNDy6HVq+OOREQk7ZQIqnLRRWAGc+fGHYmISNopEVTlyCPDiKR33aXOZSLS4ikRVGfSJFizBl54Ie5IRETSqtpEYGaDzWypma03swIz65q0bVlmwovRl78M7dvD7NlxRyIiklY1XRHMAn4ODAb+CfzDzPpG23LSHFf8OneGL30J7r0Xdu+OOxoRkbSpKREc6u6F7v6+u/8XcAVQaGanAtlx43zqVNi+HR58MO5IRETSpqZEYGbWpeKJuz8JfBn4C9A73YE1CZ/7HBx7LPzpT3FHIiKSNjUlgt8CA5IL3P1V4PNAdvyJ3KpVuCr4+9/Vp0BEWqxqE4G73+PuS6sof9vdL0tvWE3I1KmQkwN//GPckYiIpEWtzUfNLN/MHjKzl8zsVTN7zcxezURwTULPnqEF0R13wM6dcUcjItLoUulHcDdwB6F+4Fzgi9E6e3znO6HS+O67445ERKTRmdfSc9bM/uHuozMUT6rq3WopPz+foqKiOr6bw8knw0cfwYoVYfgJEZHmp8ofr1SuCK4zsz+Z2VfN7EsVS4MiMTvMzO43szfMbJWZjTSzbma22MxKonVyB7ZrzWyNma02s7EV5cXFxQwePJi8vDxmzpxJbUmtAQHDd78Lq1bB4sXpeQ8RkZikkgi+CQwDxhFuCVXcHmqIm4BCdz8eGAqsAq4Blrh7P2BJ9BwzGwhMBAZFMdxSXl4OwIwZMygoKKCkpISSkhIKCwsbGFYNJk6Ez3wGfve79L2HiEgMUkkEQ909390nu/s3o+XS+r6hmXUGTgduB3D3j939fWACUDGew2zg/OjxBGCuu+9x97XAmmXLlrFx40Z27NjByJEjMTMmTZrE/Pnz6xtW7dq2hZkz4fHHYfny9L2PiEiGpZIIlkZ/lTeWY4HNwB1m9nJ026kj0NPdNwJE6yOi/Y8E1icdX1pWVkZZWRmJRGJ/YSKRoKysrMo3LCgoID8/n/z8fDZv3lz/yGfMgE6d4D//s/6vISLSxKSSCEYDy6P7843RfLQNcBIwy91PBHYS3QaqxqcqN8ysyvoAq6YSd/r06RQVFVFUVESPHj3qFTQAhx0G3/oW/O//wptv1v91RESakFQSwTigH3AWjdN8tBQodfeK8Z3vJySGd82sF0C03pS0/1FJxydyc3NJJBKUlpYeeNHSUnJzcxsQVoq+973Qwew//iP97yUikgGpJIIvuPtbyQswo75v6O7vAOvNrH9U9HlgJbAQmByVTQYWRI8XAhPNrJ2Z9QH6jRgxgl69etGpUyeWLl2KuzNnzhwmTJhQ37BS16sXXHZZGJ76rbfS/34iImmWSiK40My+XvHEzG4BGnB/BYB/A+6ObjENA34DXA+caWYlwJnRc9x9BTCPkCwKgctbt24NwKxZs5g2bRp5eXn07duX8ePHNzCsFP3oR2Ecol//OjPvJyKSRql0KOtA+Kv8z8B4YJu7fzf9odUosx3KqnLFFXDrrfDGG9C3b+37i4jEr24dyqIOXt2ADsA04IfADuCXUXl2+8lPoE0b+MUv4o5ERKRBaro1VAwUResngcOAc5LKs1uvXuGq4K674PXX445GRKTear011ETFf2sIYOvWcFto9Gh4+OHGeU0RkfSp91hDUp3DD4drr4VHHgmT14iINENKBA01cyYkEvD978O+fXFHIyJSZ0oEDdWhQ+hcVlwMc+bEHY2ISJ2lVEcQDQndD2hfUebuT6cxrto0jTqCCvv2wahRoYNZSQkcemjjvr6ISOOoXx2BmU0DngYeA34RrX/emJE1e61awU03wTvvwK9+FXc0IiJ1ksqtoSuB4cBb7j4GOJEweqgkO+UU+OY34fe/DxPYiIg0E6kkgo/c/SMAM2vn7m8A/Ws5Jjtdf324LXT55WF6SxGRZiCVRFBqZocB84HFZrYA2JDOoJqtI46A3/wGnnxSE92LSLNRpw5lZvZZoAthmsmP0xZV7ZpWZXGy8vLQwWzNmnCLqHv39L2XiEjd1L9DmZmNNrNvuvtTwPOEWcOkKq1bQ0EBvP8+XH113NGIiNQqlVZD1wE/Aq6NinKAu9IZVLM3eHAYqnrOHHj00bijERGpUSpXBBcA5xGmlMTdNwCd0hlUi/Czn8HAgWESm+3b445GRKRaqSSCjz1UJDhANNG81KZdO7jzTti4Ea66Ku5oRESqlUoimGdmtwKHmdllwBPAbekNq4UYPjzcIrrjDliwoPb9RURiUGOrITMzIAEcT5i83oDH3H1xZsKrVtNtNVTZxx+HzmZlZfDaa9CzZ+beW0TkYFW2GmpT0xHu7mY2391PBuL+8W+e2rYNk9ecfDJcemmYt8Cq/LcQEYlFKreGlprZ8LRH0pINGgQ33ACLFoUxiUREmpBUJq9fCRwHvEVoOWSEi4Uh6Q+vWs3n1lAFdzj//NCcdOlSOOmkzMcgItmuytsRqSSC3lWVu/tbjRBUfTW/RABhasuhQ6F9+zB/QZcu8cQhItmqfj2L3f2t6Ed/N+EHeH9TUqmjww+HuXNh3bpQX6CB6USkCUilZ/F5ZlYCrAWeAtYB6i5bX6NHh1FKH3wQbrwx7mhERFKqLP4VcCrwT3fvA3weeDatUbV0V18d6gt+8AN46qm4oxGRLJdKIvjE3bcCrcyslbs/CQxLb1gtnBnMng15eXDRRbB+fdwRiUgWSyURvG9mhxKmq7zbzG4C9qY3rCzQuTPMnw8ffRSuDnbtijsiEclSqSSCCcAu4CqgEHgTODedQWWN44+He+6Bl1+GKVNUeSwisUil1dBOd9/n7nvdfba73xzdKpLG8MUvwm9/C/fdB7/4RdzRiEgWqnGICcmQ738fVq4MiSAvDy65JO6IRCSLKBE0BWZw663w1luhf0EiAZ/7XNxRiUiWSGmqynQws9Zm9rKZPRw972Zmi82sJFp3Tdr3WjNbY2arzWxsRXlxcTGDBw8mLy+PmTNnUpf5l5uctm1D34J+/eCCC+D11+OOSESyRLWJwMxeM7NXq1sa4b2vBFYlPb8GWOLu/YAl0XPMbCAwERgEjANuKS8vB2DGjBkUFBRQUlJCSUkJhYWFjRBWjA47LAxM16EDjB0brhBERNKspiuCLxJaBxVGy9ejZRFwf0Pe1MwSwDnAn5KKJwCzo8ezgfOTyue6+x53XwusWbZsGRs3bmTHjh2MHDkSM2PSpEnMnz+/IWE1Db17Q2Eh7NwZksGWLXFHJCItXLWJIGmModPc/Yfu/lq0XAOMre64FN0I/BDYl1TW0903Ru+9ETgiKj8SSO5xVVpWVkZZWRmJRGJ/YSKRoKysrMo3KygoID8/n/z8fDZv3tzA0DNgyBD461/DmETjxmnOYxFJq1TqCDqa2eiKJ2Y2Cqj3vMVm9kVgk7sXp3pIFa9RZX2AVTPhy/Tp0ykqKqKoqIgePXrUJdz4/J//Aw88AK+8EpqYqsOZiKRJKolgKvA/ZrbOzNYBtwCXNuA9TwPOi15rLnCGmd0FvGtmvQCi9aZo/1LgqKTjE7m5uSQSCUpLS/cXlpaWkpub24CwmqBzzoG774bnngu9jz/6KO6IRKQFSqVDWbG7DwWGAEPdfZi7v1TfN3T3a9094e7HECqB/+bulwALgcnRbpOBitneFwITzaydmfUB+o0YMYJevXrRqVMnli5dirszZ84cJkyYUN+wmq6vfAVuvx2eeCK0JlIyEJFGVm0/AjO7xN3vMrPvVSoHwN1/38ixXA/MM7OpwNvARdH7rDCzecBKwhhHl7du3XoRwKxZs5gyZQq7d+9m/PjxjB8/vpFDaiKmTIF9+2DqVPjyl0Mz03bt4o5KRFqIamcoM7NvufutZnZdVdvdPc7xEJrnDGUNddttMH16aE300EOhmamISOqqrEit9oogSgKtgR3u/v/SFpak7rLLoHVrmDYt1B/89a/Qsd719iIiQC11BO5eDpyXoVgkFZdeCnPmhAltzjoL3n8/7ohEpJlLZayh58zsD8D/AjsrChtSYSwNdMkl0L49fO1rYUyixx6Dnj3jjkpEmqlUEsGoaP3LpDIHzmj8cCRlF14YJre54IIwD/Ljj0OfPnFHJSLNULWVxU1cdlYWV+X550OHs5wcePRROPHEuCMSkaarbpXF+4+q1Hw0sh0odvflDQxKGmrkSPjHP8JQFKefHpqWnnlm3FGJSDOSSs/ifODbhDF/jgSmA58DbjOzH6YvNEnZgAHhyuDYY+Hss0MHNBGRFKWSCA4HTnL3q939akJi6AGcDkxJY2xSF7m58MwzcMYZoXnptdeGTmgiIrVIJREcDXyc9PwToLe77wb2pCUqqZ/OneHhh0Ons+uvh4sugg8/jDsqEWniUmk1dA+w1Mwqxv45F7jXzDoShn2QpiQnB/74Rzj++DAX8mmnwcKFYZ4DEZEqpNRqyMxOBkYTapz/4e5xN7tRq6FUFBbCxReHcYnmzdM8yCJSZauhlOYsjkYgvcndb2wCSUBSNW4cvPACdOsGX/gC3HgjNM/mwiKSRrFNXi8ZcvzxsGwZnHsuXHVV6JWsegMRSaJEkA06dw6znf361zB3LowYAStVvSMigRJBtmjVCn78Y1i8GLZuheHD4S9/iTsqEWkClAiyzRlnwMsvQ34+TJoEkyfrVpFIllMiyEa5ubBkCVx3Hdx1F5x0ErykwWRFspUSQbZq0wZ+/nP4299g1y449dTQCa28PO7IRCTDlAiy3Wc/C6++ChMmhGEpxoyBdevijkpEMkiJQEI/g3nzYPZsWL4cBg+GggL1ORDJEkoEEpiFyuPXXgvNS7/1LRg/HkpL445MRNJMiUAO1rt3aGL6hz+E0UwHDoRbb9VIpiItmBKBfFqrVnD55QeuDr79bfj852HNmrgjE5E0UCKQ6h17bLg6+NOfQt+DE06Af/93+Pjj2o8VkWZDiUBqZgZTp8KqVXD++fCzn8GwYfD003FHJiKNRIlAUtOrVxinaNEi2L07NDu95BLYuDHuyESkgZQIpG7Gj4cVK+CnP4X77oP+/eH3v9ftIpFmTIlA6u6QQ+BXv4LXX4fRo+Hqq2HIEHj00bgjE5F6UCKQ+uvXDx55JMyTvG8fnH12WDTEtUizokQgDWMG55wTrg5+9zt47rnQM/nb34Z33ok7OhFJQcYTgZkdZWZPmtkqM1thZldG5d3MbLGZlUTrrknHXGtma8xstZmNrSgvLi5m8ODB5OXlMXPmTFKZf1nSpG1b+N73Ql+Dyy+H22+HvDz4xS/ggw/ijk5EahDHFcFe4Gp3HwCcClxuZgOBa4Al7t4PWBI9J9o2ERgEjANuKY9GyJwxYwYFBQWUlJRQUlJCYWFh5j+NHKx7d7j55lChPHZsGOG0b99QtmdP3NGJSBUyngjcfaO7vxQ9/gBYBRwJTABmR7vNBs6PHk8A5rr7HndfC6xZtmwZGzduZMeOHYwcORIzY9KkScyfPz+jn0VqcNxxYXrMpUtDR7QrrwwtjG6/HT75JO7oRCRJrHUEZnYMcCLwAtDT3TdCSBbAEdFuRwLrkw4rLSsro6ysjEQisb8wkUhQVlZW5fsUFBSQn59Pfn4+mzdvbvwPItU75ZQwCc5jj0GPHjBtGgwYEKbJ3Ls37uhEhBgTgZkdCjwAfNfdd9S0axXHVlkfYPapXQGYPn06RUVFFBUV0aNHj3pGLPVmBmedBcuWwYIFcOihYaTTgQPD0Ne6QhCJVSyJwMxyCEngbnd/MCp+18x6Rdt7AZui8lLgqKTDE7m5uSQSCUqThkguLS0lNzc3/cFL/ZnBeeeFaTEfeAA6doQpU8Ito9tuUx2CSEziaDVkwO3AKnf/fdKmhcDk6PFkYEFS+UQza2dmfYB+I0aMoFevXnTq1ImlS5fi7syZM4cJEyZk7HNIA7RqBV/6UkgICxeGCubp08Mgd7/7nVoZiWRYHFcEpwHfAM4ws+XRcjZwPXCmmZUAZ0bPcfcVwDxgJVAIXN66dWsAZs2axbRp08jLy6Nv376MHz8+ho8j9WYG554LL7wQRjk9/nj4/vfDnAg//an6IYhkiDXTtvf1Djo/P5+ioqLGjEUa0wsvwG9/C/PnQ04OfOMboX/CwIFxRybSElRZkaqexdK0nHIKPPggrF4dhr+++24YNCj0SVi0SDOliaSBEoE0Tf36wS23wPr1YTKc118PQ1kMGAD//d+wo6aGZiJSF0oE0rR17w4/+QmsXQv33ANdu8LMmZCbCzNmhOk0RaRBlAikeWjbFr761dBT+cUX4aKL4M47w/DXo0fDnDlhwhwRqTMlAml+8vPhjjugtBRuuAE2bYLJk8NVwpVXwquvxh2hSLOiRCDN1+GHh+amq1fD3/4G48bBrFkwdCgMHx4ev/9+3FGKNHlKBNL8mcGYMXDvvbBhA9x0U5g68zvfCXMtT5wYWhxpbCORKikRSMvSvXuoTF6+HIqLwyB3TzwRWhwdeSRcdRUUFUHz7D8jkhbqUCYt38cfh/mU58yBv/41DHLXrx987WuhArp//7gjFMmUKjuUKRFIdnnvvdBh7Z574Mknw5XB0KHwla+EJS8v7ghF0kmJAJQIJMmGDTBvXliefz6UnXgifPnLYVC8AQPijU+k8SkRgBKBVOPtt+H+++G++0JfBQiJ4EtfggkTQpPVaua7EGlGlAhAiUBSUFYGDz0UbiE99VQY3yg3N8ylMGFCaKHUrl3cUYrUhxIBKBFIHW3ZEpqeLlgQptvcuTNMqPOFL4SWSOecE5KESPOgRABKBNIAH30UOq498gg8/HC4nQQwbFjozDZuHIwcGYbDEGmalAhAiUAaiTusWBESQmEhPPts6LB26KFwxhlw5pnhqqF/f9UtSFOiRABKBJImO3aE5qiPPhpmW/vXv0J5IhESwhlnhLqFRCLeOCXbKRGAEoFkyL/+FXo0P/EELFkC27aF8n79QkIYMwZOP131C5JpSgSgRCAx2LcvjIj65JOhjuHppw9MrNO3b0gIp58Op50WOrTpVpKkjxIBKBFIE7B3b0gMTz0VksIzz8DWrWHbEUfAqFEhKYwaBSedBO3bxxuvtCRKBKBEIE3Qvn2walWocK5Y3nwzbMvJCb2dTz01LCNGwLHH6qpB6kuJAJQIpJl4553Qw7liefFF2LUrbOvWLfR0Hj48LCefHEZWVXKQ2ikRgBKBNFN798Lrr4eE8OKLsGxZeF5eHrb36BFuI1UsQ4eG+odWGmleDqJEAEoE0oLs2gWvvAIvvRSW4uLQt6FiAp6OHcOczkOHhvXgwWHp0iXeuCVOSgSgRCAt3EcfwcqVYWKe5ctDoli+/EArJYCjjgoJ4YQTYODAsAwYEDrDSUtXZSJok+koRCSN2rc/cHuogjusXw+vvXbw8sQTYdKeCkcfHRJC//5w/PFh6d8/TPep+ocWTYlApKUzCz/yRx8dBsmrsHdv6Pi2cmW4pbRyJaxeDf/4Rxhcr8Khh4b+Df36heW440L9Q9++8JnPKEm0ALo1JCIHcw9Dcb/xRlhKSg4sa9ceqKAG6NAhNGft2xf69AnLMceEpU8f6Nw5rk8hVVMdASgRiDTIJ5/AunWhn0PlZd26g68kAA477MDVSMVy1FFhzKVEIgyxoQ5zmaQ6AhFpoJycA7eIKnMPPaTXrQtXDmvXhrqJt94KQ3Y/+2yYM7qy7t1DUujVKySG5PVnPhOWnj3D1YekhRKBiDQOs/Cj3r176PBWlQ8+gNLScOuptPTgZeNGePll2LQp9LaurHPnkBB69gz9Jnr0CENyVDzu3h0OP/zA0rGj6i9S1OwTQWFhIVdeeSXl5eVMmzaNa665Ju6QRKQ6nTqFlkkDBlS/z969IRls2ADvvhuWd9458HjTJvjnP8MVxpYtVScNCBMEHX44dO16YOnWLay7dAm3rbp0OXjp3DnE2KlTSCRZ0iGvWdcRlJeXc9xxx7F48WISiQTDhw/n3nvvZeDAgdUeqDoCkRakvDzcbtq0KdyWqrxs2xa2V6wrluR+FdUxCy2mKpaOHQ+sO3aEQw4JS4cOB9YVS/v2YV7r9u0PPG7XLiSnnJywrljatAlLTs6Bx23aQOvWYTFrzCublldHsGzZMvLy8jj22GMBmDhxIgsWLKgxEYhIC9K69YHbUXVRXh5uU23fDu+/H9bbt4eyHTvCuuLxhx+GSvCK5b33wq2sXbtg9+4D608+SctHBMKVScXy3HNhfKlG1CyvCMaNG+dbtmzhvffeY8eOHfTu3RuArVu3snPnTo4++uiD9t+8eTNbtmwBYM+ePQwbNizTIddq8+bN9OjRI+4wPkVx1Y3iqpsWFZd7uE1VeV3xuGKp/LyqpeL1Kq137drFIb1713te7OLi4sfcfVwVsXtzXNzdfd68eT516tSKpz5nzhy/4oorvCaHHHJIjdvjcvLJJ8cdQpUUV90orrpRXHXTCHFV+ZvarGtCEokE69ev3/+8tLSUXE39JyJSJ806EQwfPpySkhLWrl3Lxx9/zNy5cznvvPPiDktEpFlp1pXFbdq04Q9/+ANjx46lvLycSy+9lEGDBtV4TPe6ViplyPTp0+MOoUqKq24UV90orrpJV1zNsrIYDTEhIlIfVTYfbda3hkREpOGUCEREslyLSgSFhYX079+fvLw8rr/++k9td3fWr19PXl4eQ4YM4aWXXkr52HTGdffddzNkyBCGDBnCqFGjeOWVV/ZvO+aYYxg8eDDDhg0jv7rxW9IU19///ne6dOnCsGHDGDZsGL/85S9TPjadcd1www37YzrhhBNo3bo127ZtA9L3fV166aUcccQRnHDCCVVud3dmzpyZ8XOrtrjiOrdqiyuuc6u2uOI4twDWr1/PmDFjGDBgAIMGDeKmm2761D5pPceqa1faxJdP2bt3rx977LH+5ptv+p49e3zIkCG+YsWKg/Z55JFHvHPnzr5v3z5//vnnfcSIESkfW1+pvPazzz7r27Ztc3f3RYsW7Y/L3b13796+efPmRomlrnE9+eSTfs4559Tr2HTGlWzhwoU+ZsyY/c/T9X099dRTXlxc7IMGDapy+yOPPOLjxo3L6LmVSlxxnFupxBXHuZVKXMkydW65u2/YsMGLi4vd3X3Hjh3er1+/Kn+/GuEcq/I3Ne4f9EZbgJHAY0nPrwWurbTPrcArSc9XA71SOTadcVXavytQlvR8HdA9pu/rc8DDDf1Maf6+7gEuS/f3Fb32McDr1Wy7FfhqJs+tVOKK49xK8fvK+LlVj+8rY+dWFe+9ADgzU+dYS7o1dCSwPul5aVRWeZ/Lq9gnlWPTGVeyqcCjSc8deNzMis2sMduOpRrXSDN7xcweNbOKtrlN4vsys0OAccADScXp+r5qU13c6fyu6ipT51aqMn1upSzOc8vMjgFOBF6otClt51iz7kdQSVXNoio3M61un1SOra+UX9vMxhD+s45OKj7N3TeY2RHAYjN7w92fzlBcLwG93f1DMzsbmA/0S/HYdMZV4VzgWXffllSWru+rNnGcWynL8LmVijjOrbqI5dwys0MJyee77l55iNS0nWMt6YqgFDgq6XkC2JDiPqkcm864MLMhwJ+ACe6+taLc3TdE603AQ8CITMXl7jvc/cPo8SIgx8y6p3JsOuNKMhG4t1LM6fq+ahPHuZWSGM6tWsV0btVFxs8tM8shJIG73f3BKnZJ3zmWiftdGbqn1gb4F9AHaAu8AgyqtM85hEtjA04FlqV6bJrjOhpYA4yqVN4R6JT0+DlgXAbj+gwHOh2OAN6OvrtYv69ovy7ANqBjJr6v6DWPofp73hk/t1KMK+PnVopxZfzcSiWuGM8tA+YAN9awT9rOsRZza8jd95rZFcBjQGvgz+6+wsy+HW3/I7AIOJvwH2MX8M2ajs1gXP8XOBy4xcIEFHvdPR/oCTwUlbUB7nH3wgzGdSEww8z2AruBiR7OvLi/L4ALgMfdPXm29LR9X2Z2L6GCs7uZlQLXATlJMWX83EoxroyfWynGlfFzK8W4IMPnVuQ04BvAa2a2PCr7MSGRp/0ca65DTIiISCNpSXUEIiJSD0oEIiJZTolARCTLKRGIiGQ5JQIRkSbCzP5sZpvM7PUU9j3dzF4ys71mdmGlbZPNrCRaJtf2WkoE0uKZWbmZLU9arqll/x838P3ON7OB1Wy7s/J/2jq87jAze97MVpjZq2Z2cUPilCbpTsLQFql4G5hCGBNpPzPrRmgWewqhj8Z1Zta1phdSIpBssNvdhyUttY3T26BEAJwPVJkIGmgXMMndBxF+LG40s8PS8D4SEw9DViQPa4GZ9TWzwmiMo2fM7Pho33Xu/iqwr9LLjAUWu/s2d38PWEwtyUWJQLKSmXUxs9Vm1j96fq+ZXWZm1wMdoiuHu6Ntl5jZsqjsVjNrHZV/aGa/jgZOW2pmPc1sFHAecEO0f98aYvhVdIXQyszWmdlvor/4i8zsJDN7zMzeTOpM9093L4kebwA2AT3S+kVJU1AA/Ju7nwx8H7illv3rPAidEoFkg4of9orlYnffDlwB3GlmE4Gu7n6bu1/DgSuIr5vZAOBiwoBjw4By4OvR63YElrr7UOBpwpDFzwELgR9Er/FmVQGZ2X8CRwDfdPeKv+jWu/tI4BnCLYILCUMJ/LKK40cQhhOo8vWlZYgGoRsF3Bf1OL6VMPR0jYdVUVZjz+EWM8SESA12Rz/iB3H3xWZ2EfA/wNBqjv08cDLwYjS8QAfCX+IAHwMPR4+LgTNTjOdnwAvuXnko44XR+jXgUHf/APjAzD4ys8Pc/X0AM+sF/AWYnJREpGVqBbxf1flbg1LCMBoVEsDfa3sTkaxkZq2AAYSxbrpVtxswO6l+ob+7/zza9okfGKOlnNT/sHoRODmq1Eu2J1rvS3pc8bxNFHNn4BHgp+6+NMX3k2bKw1DUa6M/WLCguj9aKjwGnGVmXaNK4rOismopEUg2uwpYBXwV+HM0DDDAJ0mPlwAXWhiDHjPrZma9a3ndD4BONWwvBK4HHjGzmvY7iJm1JQx/PMfd70v1OGk+okHxngf6m1mpmU0l3IqcamavACuACdG+w6OB8y4CbjWzFQAe5lD4FeEPjheBX/rB8yp8im4NSTbokDSiI4Qf4j8D04AR7v6BmT0N/JTQ7K4AeNXMXorqCX5KmJmqFfAJYZa7t2p4v7nAbWY2E7iwqnoCd78vSgILLUzMkoqvAKcDh5vZlKhsirsvr/YIaVbc/avVbPpUqx93f5Fw26eq1/kz4RxPiUYfFRHJcro1JCKS5ZQIRESynBKBiEiWUyIQEclySgQiIllOiUBEJMspEYiIZLn/D2W0N+F6cKxXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid_size_km = np.linspace(100,1,100) # 100km to 1 km\n",
    "width = np.linspace(1000,100000,100) # 100000km to 1 km \n",
    "height = np.linspace(1000,200000,100)# 200000km to 1 km\n",
    "areaExtent = (width)*(height)  #153531695409.90253 m2 the geometry above spain\n",
    "areaGrid = (grid_size_km) ** 2\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "def calculateSquareOverExtent(geometry):\n",
    "    \"\"\"Calculates the number of squares in grid covered by an extent and a grid_size \n",
    "    Parameters:\n",
    "    geometry (List(float)): [x_min, x_max, y_min, y_max] in epsg 4326\n",
    "    grid_size_km (float):  grid cell size in km\n",
    "\n",
    "    Returns:\n",
    "    int: Returning an aproach number of grids cells\n",
    "    \n",
    "    \"\"\"\n",
    "    wgs84 = pyproj.CRS('EPSG:4326')\n",
    "    webMercator = pyproj.CRS('EPSG:3857')\n",
    "    project = pyproj.Transformer.from_crs(wgs84, webMercator, always_xy=True).transform\n",
    "    utm_extent = transform(project, geometry).bounds\n",
    "    width = utm_extent[2] - utm_extent[0]\n",
    "    height = utm_extent[3] - utm_extent[1]\n",
    "    areaExtent = width*height\n",
    "    min_grid = calcMinGrid(areaExtent)\n",
    "    max_grid = calcMaxGrid(areaExtent)\n",
    "    grid_size_m = round(min_grid/1000) * 1000\n",
    "    areaGrid = grid_size_m ** 2\n",
    "    print(f'geometry extent: {round(areaExtent/1000,2)} km2')\n",
    "    print(f'min grid size: {round(min_grid/1000)}')\n",
    "    print(f'max grid size: {round(max_grid/1000)}')\n",
    "    return round(areaExtent/areaGrid)\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.xlabel(\"Extent km2\")\n",
    "plt.ylabel(\"grid area km2\")\n",
    "\n",
    "# plot the function\n",
    "plt.plot( areaExtent, areaGrid,'r')\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "7f7edf57",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAr3UlEQVR4nO3deXxU5fX48c8RFBW1oiwGAoIEI2Aj1oiCG8oiRQk/KkaoC4o21n2rFq3WpVpiq34FjGJkCwWJBZVQWVJEUdzAiKKCShQQEgKERQWUsOT8/ngGjGmWSzIzd5bzfr3ySmbmzr3nKhyenHme84iqYowxJnYc5HcAxhhjgssSuzHGxBhL7MYYE2MssRtjTIyxxG6MMTGmod8BAHWeltO3b1/mzp0bzFiMMSZaSHUvRPWIfdOmTX6HYIwxESeqE7sxxpj/FbLELiLjRWSjiHxe6flbROQrEVkmIv/Y9/yIESNISkoiOTmZ/Pz8UIVljDExL5Q19onAM8CkfU+IyPnAACBFVctEpDlw9/Lly8nNzWXZsmWsW7eOXr16sWLFCho0aBDC8IwxJjaFbMSuqm8DWyo9fQOQqaplgWM2AuTl5TF48GAaNWpEu3btSEpKYvHixaEKzRhjYlq4a+wnAueIyCIReUtETgcoLi6mdevW+w9KTEykuLi4yhNkZ2eTmppKamoqpaWlYQnaGGOiSbgTe0OgCXAmcDfwb1WlqkZkIlXP5MnIyKCgoICCggKaNWsWyliNMSYqhTuxFwGvqLMYKN+0aROJiYmsXbv254OKimjZsmWYQzPGmNgQ7sQ+A7gAQEROBA5p2rQpaWlp5ObmUlZWxqpVqygsLKRr165hDs0YY8Lk00/h2mvhvfdCcvqQzYoRkalAD6CpiBQBDwLjgfGBKZC7gKEiMr9z586kp6fTqVMnGjZsSFZWls2IMcbElr17YeZMGDkS3noLDj8cunWD7t2DfimJgI026hxAamoqBQUFwYzFGGOC67vvYNw4eOYZWL0a2rSBm2+G666DJk3qc+ZqWwpEQq8YY4yJPV99BaNGQU4O7NgB55wDTzwBAwZAw9CmXkvsxhgTLOXlkJ/vEvrcuXDIIfD738Ott8Kpp4YtDEvsxhhTX9u3u5H56NFupH7ccfDII3D99dC8edjDscRujDF1tWqVq52PGwfffw9du8KUKTBokBut+8QSuzHGHAhVWLDAzW6ZORMaNHCJ/Lbb4Mwz/Y4OsMRujDHe/PQTvPiiq59/+ikceyzcey/ceCO0auV3dL9gid0YY2pSVATPPQfPPw+bN0NKiiu9DBkChx3md3RVssRujDGVqcIHH7hyy8svu8VFAwa4cst550E1vawihSV2Y4zZZ9cumDbNJfQPP4Rf/cpNVbz5ZmjXzu/oPLPEbowxGzfCmDGu5LJ+PSQnQ1YWXHUVHHGE39EdMEvsxpj49fHHbnQ+daobrfft68otffrAQdG7JbQldmNMfNmzB/LyXEJfuBAaN3Z9W265BU46ye/ogsISuzEmPmzZAmPHuhLLmjXQti08+SQMGwZHH+13dEFlid0YE9uWL3dzzydNcnPRe/Rwo/X+/d3iohhkid0YE3vKy2HOHJfA582DRo3giivcDJeUFL+jCzlL7MaY2LFtG0yY4Jpxff01tGwJjz0GGRnQtKnf0YVNyD72FZHxIrIxsFtS5df+JCIqIvv/S48YMYKkpCSSk5PJz88PVVjGmFj0zTdw++1uaf9tt7kkPnWq29jivvviKqlDaEfsE4FngEkVnxSR1kBvYM2+55YvX05ubi7Lli1j3bp19OrVixUrVtj2eMaY6qnC/Pmu3DJrltu8Ij3dlVvifM/kkI3YVfVtYEsVL/0fcA8VtsTLy8tj8ODBNGrUiHbt2pGUlMTixYtDFZoxJpr9+CNkZ8Ovfw29e8OiRXD//fDttzB5ctwndQhzjV1E0oBiVV0qFXotFBcXc2aFdpeJiYkUFxdXeY7s7Gyys7MBKC0tDWm8xpgIsnatm6r4wgtu6mKXLq6ePngwHHqo39FFlLAldhE5HPgL0Kfya1VtqC3VNNnJyMggIyMDcJtZG2NimCq8+64rt7z6qns8cKCro599dsQ34/JLOEfs7YF2wL7ReiKwZP369SQmJrJ27dr9BxYVFdGyZcswhmaMiShlZfDSSy6hL1niFhDdeSfcdBMcf7zf0UW8sDVDUNXPVLW5qrZV1bZAEfCb4447jrS0NHJzcykrK2PVqlUUFhbS1epkxsSf9evhoYdc8h461C0oGjPG9UT/xz8sqXsUshG7iEwFegBNRaQIeFBVx1V1bOfOnUlPT6dTp040bNiQrKwsmxFjTDwpKHCj85degt274aKLXLmlVy8rt9SBVFXfDrM6B5CamkpBQUEwYzHGhMvu3a5uPnIkvPeea497zTWuGVeHDn5HFw2q/RfPVp4aY8Jr82Y3syUry5VY2reHp5+Gq692G1uYerPEbowJj88/d6PzyZNh507o2ROefRb69YvZZlx+scRujAmdvXvdqtCRI+GNN9x88yuvdKtDTz7Z7+hiliV2Y0zwff89jB8PzzwDK1dCYiJkZroNLY491u/oYp4ldmNM8KxY4TorTpwI27fDWWe5hD5woOvlYsLC/ksbY+pH1fU8HzkSZs+Ggw92y/xvuw1OO83v6OKSJXZjTN3s2OF2JRo9Gr74Alq0cIuLrr8ejjvO7+jimiV2Y8yB+fZbVzsfOxa++86NyidNci1zGzXyOzqDJXZjjBeqsHChK7fMmOFWg15yiZvd0r27rQ6NMJbYjTHV27nT7UQ0ahR88gkccwzccw/ceCO0bu13dKYaltiNMf9r3Tp47jl4/nkoLYXOnd3mFpdfDocf7nd0phaW2I0xP1u0yJVbpk1zi4suvtjNbrngAiu3RBFL7MbEu927Yfp0l9AXLYKjjoKbb3Zf7dv7HZ2pA0vsxsSr0lJXXnn2WVd66dDB1dKvvhqOPNLv6Ew91JjYReRQ4GLgHKAl8BPwOTBLVZeFPjxjTNAtXepG5y++6HYq6tPHdVvs2xcOCtveOyaEqk3sIvIQ0B9YACwCNgKHAicCmYGkf5eqfhr6MI0x9bJ3L8yc6RL6W2+5D0D39T7v1Mnv6EyQ1TRi/1BVH6rmtadEpDnQpro3i8h43Gh/o6qeHHjun7h/LHYB3wDX7NvoY8SIEYwbN44GDRowatQoLrzwwgO+GWNMJd99B+PGuQVFq1dDmzZui7nrroMmTfyOzoRItYldVWfV9EZV3YgbxVdnIvAMMKnCc/OAe1V1j4g8DtwLsHz5cnJzc1m2bBnr1q2jV69erFixwrbHM6auvvzS1ctzcuDHH+Hcc+HJJyEtzZpxxYFa/w+LSDPgz0AnXCkGAFW9oKb3qerbItK20nP/rfDwA2AQQF5eHoMHD6ZRo0a0a9eOpKQkFi9eTLdu3bzfiTHxrrwc8vNduSU/Hw45BH7/e7c69NRT/Y7OhJGXT0qmAF8A7YCHgdXAh0G49jBgDkBxcTGtK6xiS0xMpLi4uMo3ZWdnk5qaSmpqKqWlpUEIw5got327K7V07Oh2I1q6FB55BNauhQkTLKnHIS+J/VhVHQfsVtW3VHUYcGZ9LioifwH24P7RoKoNtaWaxRAZGRkUFBRQUFBAs2bN6hOGMdFt5Uq4805o1cp9CHr00TBlimvS9cAD0Ly53xEan3gptu0OfC8RkYuAdUBiXS8oIkNxH6r21EBGT0xMZO3atfuPKSoqomXLlnW9hDGxSxUWLHDllpkz3V6hgwa51aFn1mu8ZWKIlxH7oyLyK+Au4E/AWOCOulxMRPri6vVpqvrjvufT0tLIzc2lrKyMVatWUVhYSNeuXetyCWNi008/uTa5p5zilve/+y7cd5+b6TJ1qiV18wu1jthV9bXAj98D53s9sYhMBXoATUWkCHgQNwumETAvUGr5QFXp3Lkz6enpdOrUiYYNG5KVlWUzYowBKCpyK0Ozs2HzZkhJcdMXhwyBww7zOzoToaSq+vYvDhA5EXgOaKGqJ4tICm7E/WiQYqg5gBqkpqZSUFAQpDCMiRCq8MEHrtwyfbp7nJbmyi3nnWfNuMw+1f5B8FKKeQE30t4NEFhpOjg4cRlj9tu1CyZPhq5d3eYVc+fC7bfD11/Dq69Cjx6W1I0nXj48PVxVF1eapbInRPEYE382bHB9z597Dtavh+RkyMqCq66CI47wOzoThbwk9k0i0p5AyUREBgElIY3KmHjw8ceu3DJ1qhut9+3ryi19+lgzLlMvXhL7TUA2cJKIFAOrgMtDGpUxsWrPHrdn6MiR8M470Lgx/OEPbh56crLf0ZkYUVvb3gbADaraS0QaAwep6rbwhGZMDNmyxU1XzMqCNWugbVvXu2XYMLewyJggqjGxq+peETkt8POO8IRkTAxZtsw14/rXv9xc9B493Gi9f3+3uMiYEPBSivlYRGYC04D9yV1VXwlZVMZEs/JymD3bJfDXX4dGjeCKK1wzrpQUv6MzccBLYj8G2AxU7OaogCV2Yyr64QeYOBFGj3ZTFFu1gsceg4wMaNrU7+hMHPGy8vSacARiTNT6+muXzCdMgG3boFs3ePRR+N3v4OCD/Y7OxCEv/dgPBa4FOvPLfuzDQhiXMZFNFebPd+WWWbPc5hXp6W664umn+x2diXNeJsv+CzgOuBB4C9fZ0WbGmPj044+ub8uvfw29e8OiRXD//a5V7uTJltRNRPBSY09S1UtFZICq5ojIi0B+qAMzJqKsWeOmKr7wAmzd6javmDgRLrsMDj201rcbE04H0o/9OxE5GVgPtA1ZRMZEClXXHnfkSNerRRUGDnTllrPPtr4tJmJ5SezZItIEuB+YCRwBPBDSqIzxU1kZ5Oa6+edLlrgFRHfeCTfdBMcf73d0xtTKy6yYsYEf3wZOCG04xvho/XrXiGvMGNi40e0hOmaMm4PeuLHf0RnjmZcRuzGxraDAlVteegl274aLLnLlll69rNxiolLIWsiJyHgR2Sgin1d47hgRmScihYHvTfa9NmLECJKSkkhOTiY/3z6bNSG2e7dL5N27u5ksM2bADTfAihXw2mtuxosldROlQtkbdCLQt9Jzw4H5qtoBmB94zPLly8nNzWXZsmXMnTuXG2+8kb1794YwNBO3Nm+GESPghBNg8GBXcnn6aSgudqP2Dh38jtCYequ2FCMiv6vpjbX1ilHVt0WkbaWnB+D2QQXIARYA5OXlMXjwYBo1akS7du1ISkpi8eLFdOvWrZbwjfHos8/ch6GTJ8POndCzp9tLtF8/a8ZlYk5NNfb+ge/Nge7AG4HH5+MScl16xbRQ1RIAVS0RkeYAxcXFnFlhl/XExESKi4urPEF2djbZ2dkAlJaW1iEEEzf27nVllZEj4c033ebPV13lep+ffLLf0RkTMtUm9n09YkTkNaDTvoQsIglAVjCDqGpDbammvpmRkUFGRgbgNrM25n98/z2MHw/PPAMrV0Lr1pCZCdddB8ce63d0xoScl1kxbfcl9YANwIl1vN4GEUkIjNYTgI3A0YmJiaxdu3b/QUVFRbRs2bKOlzBxa8UK14xr4kTYvh3OOssl9IEDXS8XY+KElw9PF4hIvohcLSJDgVnAm3W83kxgaODnoUAeQFpaGrm5uZSVlbFq1SoKCwvp2rVrHS9h4ooq5Oe7WnlystsUeuBAN4XxnXfg0kstqZu442WB0s2BD1LPCTyVraqv1vY+EZmK+6C0qYgUAQ8CmcC/ReRaYA1wKXB3586dSU9Pp1OnTjRs2JCsrCwa2AdapiY7dsCkSe4D0S+/hBYt4KGH4Prr4bjj/I7OGF9JVfXtMKtzAKmpqRQUFAQzFhPpVq92tfNx4+C77+C009xiovR0t1ORMfGj2oUWNU13fEdVzxaRbfwy+QqgqnpUEAM0pnqq8PbbbnZLXp5bOHTJJS6hd+tmC4mMqaSmWTFnB74fGb5wjKlg506YOtUl9KVL4Zhj4J574MYb3UwXY0yVaqyxi8hBwKeqapN+TfisW+cWDz3/PGza5OacZ2fD5ZfD4Yf7HZ0xEa/GxK6q5SKyVETaqOqacAVl4tSiRW50Pm2aW1zUv78rt5x/vpVbjDkAXuaBJQDLRGQxsGPfk6qaFrKoTPzYtQteftkl9EWL4Kij4Oab3Vf79n5HZ0xU8pLYHw55FCb+lJa6Usuzz0JJiWu+NWoUXH01HGkf6xhTH17msb8VjkBMnFi61I3OX3zR7VTUpw+MHQt9+8JBoWw2akz8qDWxVzHdEeB7oAC4S1VXhiIwE0P27nXTFEeNgrfech+AXnMN3Hqr26XIGBNUXkoxTwHrgBdxc9gHA8cBXwHj+bkNrzG/tHWrW0j0zDPw7bduv9B//hOuvRaaNKn9/caYOvGS2Puq6hkVHmeLyAeq+oiI3BeqwEwU+/JLNzrPyYEff4Rzz4WnnoK0NOvbYkwYePlbVi4i6cD0wONBFV7zvR+BiRDl5TB3rquf//e/bnn/73/vyi1duvgdnTFxxUtivxwYCTyLS+QfAFeIyGHAzSGMzUSDbdvcyHz0aNc2NyEB/vY3yMiA5s39js6YuORlVsxKft5NqbJ3ghuOiRorV/7cjOuHH6BrV5gyBQYNgkMO8Ts6Y+KaFTyNd6pui7mRI+E//3F7hV56qSu3VNja0BjjL0vspnY//eRG46NGuU2hmzaF++6DG26AVq38js4YU4kldlO9oiK3MjQ7GzZvhpQUV3oZMsRtDG2MiUheFig1Ai4B2lY8XlUfCV1Yxjeq8P77rtzy8svucVqaa8Z13nnWjMuYKOBlDXceMADYg2sCtu+rzkTkDhFZJiKfDxkyhJ07d7JlyxZ69+5Nhw4d6N27N1u3bq3PJcyB2rULJk92H4KedZbbR/T22+Hrr+HVV6FHD0vqxkSJWrfGE5HPg9mPXURa4WbTdFLVn9LT07Vfv34sX76cY445huHDh5OZmcnWrVt5/PHHazyXbY0XBBs2wJgx7mv9erch9K23wlVXwRFH+B2dMaZ61Y60vIzY3xORXwcxGHAlncNEpOGPP/5Iy5YtycvLY+jQoQAMHTqUGTNmBPmS5heWLIGhQ6FNG7cJ9KmnugVGy5e7HYosqRsTtbx8eHo2cLWIrALK+HnP05S6XFBVi0XkCWAN8NOvfvUr+vTpw4YNG0hISAAgISGBjRs3Vvn+7OxssrOzASgtLa1LCPFrzx5XVhk1Ct55Bxo3hj/8AW65xY3UjTExwUti/20wLygiTXA1+3bAdzt27Ng1efJkz+/PyMggIyMDcKUY49Ebb7iOimvWQLt28OSTMGwYHH2035EZY4Ks2sQuIkep6g/AtiBfsxewSlVLASZNmsR7771HixYtKCkpISEhgZKSEprbcvTgycmB666DE0+EGTPg4ovd4iJjTEyqqcb+YuD7R7je6x9V+KrPJ5ZrgDNF5HARkfnz59OxY0fS0tLIyckBICcnhwEDBtTjEgZwUxX/+le3K1GPHvDeezBggCV1Y2JctSN2Vb048L1dMC+oqotEZDqwBNhTXl5ORkYG27dvJz09nXHjxtGmTRumTZsWzMvGn7IyN0qfPNmVXMaMgYMP9jsqY0wY1DrdEUBEUvjfBUqvBCmGOrf+temO1di6FQYOdLsVPfqoW/5vc9CNiTXV/qX2svJ0PJACLAPKA08rEKzEboJp5Uro1w9WrXL7ig4Z4ndExpgw8zIr5kxV7RTySEz9LVoE/fu7aY2vvw7nnON3RMYYH3hZoPS+iFhij3SvvOI+ID3ySNfrxZK6MXHLy4g9B5fc1xOEBUomyFTh6afhrrvgjDNg5kxo1szvqIwxPvKS2McDVwKf8XON3USCPXtco66sLLdz0aRJ1k7XGOMpsa9R1Zkhj8QcmO3b3Qejr70Gd98NmZlwkJfKmjEm1nlJ7F+KyIvAf3ClGCCo0x3NgVq3zq0eXbrUbYRxww1+R2SMiSBeEvthuITep8JzNt3RL599BhddBFu2uH1H+/XzOyJjTISpNbGr6jXhCMR4MG+eq6UfcQQsXOha7RpjTCXVFmVF5H4ROaaG1y8QkYtDE5b5H+PHu9F527ZuvroldWNMNWoasX8G/EdEduL6upQChwIdgC7A68DfQx1g3FOFBx6Axx6DCy+Ef/8bjjrK76iMMRGspiZgeUCeiHQAzgISgB+AyUCGqv4UnhDjWFmZ66E+darbECMryxp5GWNq5aXGXggUhiEWU9Hmza6R18KFbirjPfdYIy9jjCdeZsWYcPvmG1dP//ZbyM2Fyy7zOyJjTBSxxB5p3n8f0tKgvNw18jr7bL8jMsZEGVuqGEmmT4fzz3f7kH7wgSV1Y0yd1LTn6Whq2ARDVW+t60VF5GhgLHDySSedxPjx40lOTuayyy5j9erVtG3bln//+980adKkrpeILqrwxBOujt69O+TlQdOmfkdljIlSNY3YK+9zWvmrPkYCc1X1pKVLl9KxY0cyMzPp2bMnhYWF9OzZk8zMzHpeIkrs2QM33uiSeno6zJ9vSd0YUy+etsYDEJHGqrqj3hcUOQpYCpyg7uIKkJyczIIFC0hISKCkpIQePXrw1Vdf1XiuqN8ab9s2GDwYZs+GP/8Z/v53a+RljPGq2mlytWYREekmIsuBLwKPTxGRZ+sRzAm4xU4TROTj6667jh07drBhwwYSEhIASEhIYOPGjVW+OTs7m9TUVFJTUyktLa1HGD4rLoZzz4X8fMjOtu6Mxpig8ZJJngYuBDYDqOpS4Nx6XLMh8BvgOVU9tXHjxgdUdsnIyKCgoICCggKaReuGEp9+6jbF+PprmDXLLT4yxpgg8TREVNW1lZ7aW49rFgFFqroIYNCgQSxZsoQWLVpQUlICQElJCc2bN6/HJSJYfv7Ps13eece1CTDGmCDyktjXikh3QEXkEBH5E4GyTF2o6vrAOZMB5s+fT6dOnUhLSyMnJweAnJwcBgwYUNdLRK7sbNdy94QTXCOvU07xOyJjTAzyskDpj7hZLK1wo+3/AjfV87q3AFNE5JABAwYwYcIEysvLSU9PZ9y4cbRp04Zp06bV8xIRpLwc/vIXV0fv29c18jrySL+jMsbEqFpnxYhIc1XdWOm5ZFWtecqKd96m5VQhKmbF7NwJV18NL70Ef/wjjB4NDW3BrzGm3uo+KwZYKCLp+88kchfwajCiinmbNkGvXi6p/+Mfbhs7S+rGmBDzkmV6ANkicinQAldf7xrKoGJCYaFr5LV2rSu9XHqp3xEZY+JErSN2VS0B5gLdgLbAJFXdHuK4otu770K3bvDdd/DGG5bUjTFh5WWB0jzgDOBkoB/wfyLyRKgDi1ovvQQ9e8Ixx7hGXt27+x2RMSbOeKmxZ6nqVar6nap+DnQHvg9xXNFHFR5/3LUIOP101363fXu/ozLGxCEvpZgZlR7vUdW/hSyiaLR7N1x/PQwf7hL7vHlw7LF+R2WMiVPVJnYReSfwfZuI/BD4vu/rh/CFGOF++AH694cXXoD77oMpU+DQQ/2OyhgTx2razPrswHdbSVOdoiK3knTZMhg7Fq691u+IjDHG29Z4IvIb4GzcYqJ3VPXjkEYVDT75xCX1bdtc290+ffyOyBhjAG+zYv4K5ADHAk2BiSJyf6gDi2hz5sA550CDBm5qoyV1Y0wE8TJiHwKcqqo7AUQkE1gCPBrKwCLWmDFw882QkgKvvQYtW/odkTHG/IKX6Y6rgYqfBjYCvglJNJGsvNxtX3fDDa6R19tvW1I3xkQkLyP2MmBZYKGSAr2Bd0RkFNRvU+uo8dNPcNVVMH2625905Ejr+WKMiVhestOr/LLp14LQhBKhSkthwAC3ivTJJ+GOO0CqbapmjDG+qzaxi0g2MAd4RVW3hS+kCPLVV66R17p1MG0aXHKJ3xEZY0ytahqxjwf6AneKyC7cBhtzA3uexr6FC91IvWFDePNNOPNMvyMyxhhPqv3wVFU/UNWHVPUcIB1YA9wlIp+IyPiKPdrrQkQaiMjHF198MQBbtmyhd+/edOjQgd69e7N169b6nL5+XnzR9VFv3tyVYCypG2OiiNfNrDer6tRAM7AuQBbQoZ7Xvo0Ke6dmZmbSs2dPCgsL6dmzJ5mZmfU8fR2owt//Dpdf7pL5e++5/UmNMSaK1PrhqYgcDVyF68W+//j6zIYRkUTgIuAx3Dx58vLyWLBgAQBDhw6lR48ePP7443W9xIHbvdtNZRw3ziX2ceOgUaPwXd8YY4LEy6yY2cAHwGdAeZCu+zRwD7C/D82GDRtISEgAICEhgY0bN1b5xuzsbLKzswEoLS0NTjTbtrkPRufNgwcegIcftpkvxpio5SWxH6qqdwbrgiJyMbBRVT8SkR4H+v6MjAwyMjIAt5l1UPz5z26no/Hj4ZprgnNOY4zxiZca+79E5A8ikiAix+z7qsc1zwLSRGQ1kPvGG29wxRVX0KJFC0pKSgAoKSmhefPm9bjEAfj0U3j+ebjpJkvqxpiY4CWx7wL+CbwPfBT4KqjrBVX1XlVNVNW2wOALLriAyZMnk5aWRk5ODgA5OTkMGDCgrpc4kGDg9tuhSRN46KHQX88YY8LASynmTiBJVTeFMpDhw4eTnp7OuHHjaNOmDdOmTQvl5ZxXX3Vz1J991iV3Y4yJAaKqNR8gMhMYrKo/hiiGmgOoQWpqKgUFdfzlYedO6NgRjjgCPv7Yer8YY6JNtTM8vGSzvcAnIvImriEYEAPNv556Clavhtdft6RujIkpXjLajMBX7Fi3zi1EGjgQevb0OxpjjAkqL4n9c1X9qOITItI/RPGEx733ugVJTzzhdyTGGBN0XmbFvCAiv973QESGANG7Nd6iRTBpEtx1l7ULMMbEJC8j9kHAdBG5HLeh9VVAdG7yWV4Ot90Gxx3nRu3GGBODak3sqrpSRAbj6uxrgT6q+lOoAwuJKVPciH3iRDjyyFoPN8aYaFTtdEcR+YxfTkVsDnxPYGaMqqYEKYbwTHfcvh1OPBESE10r3oM8NbY0xphIVafpjheHIBD/jBgBJSXw8suW1I0xMa3axK6q34YzkJBatcrtV3rFFdCtm9/RGGNMSMXH0PXuu6FBA/Bj8w5jjAmz2E/sb77pyi/33gutWvkdjTHGhFxsJ/Y9e1z3xuOPd/PWjTEmDsR2k5SxY12/9WnT4LDD/I7GGGPCInZH7Fu3wv33w3nnuW3vjDEmTsRuYn/4YdiyBZ5+2vYvNcbEldhM7F98AVlZ8Ic/QJcufkdjjDFhFfbELiKtReRNEflCRJaNHDkSgC1bttC7d286dOhA79692bp1a90uoAp33AGNG8OjjwYxcmOMiQ5+jNj3AHepakfgzKysLJYvX05mZiY9e/aksLCQnj17klnXOeezZ0N+Pjz4IDRrFsy4jTEmKoQ9satqiaouCfy8rWPHjhQXF5OXl8fQoUMBGDp0KDNmzDjwk+/a5Ubryclw003BDNsYY6KGr9MdRaRt69atOeOMM9iwYQMJCQkAJCQksHHjxgM/YVYWFBbCrFlwyCFBjtYYY6KDb4ldRI4AXn766ac56qijPL8vOzub7OxsAEpLS3/54oQJcNZZ0K9fECM1xpjo4susGBE5GHgZmPK73/0OgBYtWlBSUgJASUkJzZs3r/K9GRkZFBQUUFBQQLOKNfS1a+Gzz+D//b/QBm+MMRHOj1kxAowDvlDVp/Y9n5aWRk5ODgA5OTkMGDDgwE48Z477bqN1Y0yc82PEfhZwJXCBiHzSpUsXZs+ezfDhw5k3bx4dOnRg3rx5DB8+/MDOOnu26wnTsWMoYjbGmKgR9hq7qr7DL3f+2L+D0vz58+t20rIyeP11uPJKW2VqjIl7sbHydOFC2LHDyjDGGEOsJPbZs930xgsu8DsSY4zxXWwk9jlzoEcP10bAGGPiXPQn9pUr4csvrQxjjDEB0Z/YbZqjMcb8QvQn9tmzISkJOnTwOxJjjIkI0Z3Yy8vhjTdstG6MMRVEd2LfsQN27oQ+ffyOxBhjIkZ0J/adO933k0/2Nw5jjIkg0Z3Yy8qgUSNo3drvSIwxJmJEd2LfuRPat4eDovs2jDEmmKI7I5aV2WwYY4ypJHoTe3m5JXZjjKlC9Cb24mJQdXPYjTHG7Be9iX3fnqjHHedvHMYYE2GiN7Fv2eK+H3OMv3EYY0yEscRujDExJuIS+9y5c0lOTiYpKYnMzMzqD9yxw30/4ojwBGaMMVEiohL73r17uemmm5gzZw7Lly9n6tSpLF++vOqD9606PfTQ8AVojDFRIKIS++LFi0lKSuKEE07gkEMOYfDgweTl5VV98N697nuDBuEL0BhjooCoau1HhVDfvn1106ZNAGzdupUffviB448/HoDNmzezY8cO2rRps//40tJSNm3axOF797Jt1y5O7tIlblaelpaW0qxZM7/DCAu719gVT/cbynv96KOP8lW1b1Wv+Z7Ygf0BTJs2jfz8fMaOHQvAv/71LxYvXszo0aOrfGPjxo3Zsa/WHgdSU1MpKCjwO4ywsHuNXfF0vyG+V6nuhYga6iYmJrJ27dr9j4uKimjZsqWPERljTPSJqMR++umnU1hYyKpVq9i1axe5ubmkpaX5HZYxxkSVhn4HUFHDhg155plnuPDCC9m7dy/Dhg2jc+fO1R7ftGnTMEbnv4yMDL9DCBu719gVT/fr171GVI39QMVTrc4YYyqJjhq7McaY+rPEbowxMSYqEnttbQZUlVtvvZWkpCRSUlJYsmSJD1EGR233OmXKFFJSUkhJSaF79+4sXbrUhyiDx2sLiQ8//JAGDRowffr0MEYXXF7udcGCBXTp0oXOnTtz3nnnhTnC4KntXr///nv69+/PKaecQufOnZkwYYIPUQbHsGHDaN68OSdXs/eyL/lJVf3+qtGePXv0hBNO0G+++UbLyso0JSVFly1bpqqqp512mqqqzpo1S/v27avl5eX6/vvva9euXWs7bUSq6V73effdd3XLli2qqjp79uyovVdVb/e777jzzz9ff/vb3+q0adN8iLT+vNzr1q1btWPHjvrtt9+qquqGDRv8CLXevNzrY489pvfcc4+qqm7cuFGbNGmiZWVlfoRbb2+99ZZ+9NFH2rlz5ypfD2F+qjav+p3Ua/0CugH5FR7fC9wb+Hlu4PvzwJAKx3wFJPgdezDvtZrjmwDFfscd6vsFbgduAiYCg/yOO1T3CtwIPOp3rGG613uBZ3EfALYDvgYO8jv2etxzW+Dzal4Le36KhlJMK2BthcdFgefQn5fTVntMlDnQ+7gWmBPSiEKr1vsVkVbAQGBMGOMKBS//b08EmojIAhH5SESuClt0weXlXp8BOgLrgM+A21S1PDzhhV3Y81NEzWOvRlVTeipPkfRyTDTwfB8icj4usZ8d0ohCy8v9Pg38WVX3ilQ7uysaeLnXhsBpQE/gMOB9EflAVVeEOrgg83KvFwKfABcA7YF5IrJQVX8IcWx+CHt+iobEXgS0rvA4Efev/IEeEw083YeIpABjgd+q6uYwxRYKXu43FcgNJPWmQD8R2aOqM8ISYfB4/XO8SVV3ADtE5G3gFCDaEruXe70GyFRXm/haRFYBJwGLwxNiWIU9P0VDKeZDoIOItBORQ4DBwMxKx8wErhLnTOB7VS0Jd6BBUOu9ikgb4BXgyigcyVVW6/2qajtVbauqbYHpwI1RmNTB25/jPOAcEWkoIocDZwBfhDnOYPByr2twv5kgIi2AZGBlWKMMn7Dnp4gfsavqHhG5GcgHGgDjVXWZiPwx8PoYYDbQD/cBzI+40UDU8XivfwWOBZ4NjGL3qGqqXzHXh8f7jQle7lVVvxCRucCnQDkwVlU/9y/quvH4//VvwEQR+QxXqvizqm7yLeh6EJGpQA+gqYgUAQ8CB4N/+SkSWgoYY4wJomgoxRhjjDkAltiNMSbGWGI3xpgYY4ndGGNijCV2Y4wJIxEZLyIbRaTWGU8i8n8i8knga4WIfOflGpbYTUQSkaNF5MYKj1uKyPTAzz1E5LUQXruLiPQLwXkfEZFewT6viToTgb61HQSgqneoahdV7QKMxq1hqZUldhOpjsY1xQJAVdep6qAwXbsLbt5xUKnqX1X19WCf10QXVX0b2FLxORFpLyJzAz2CForISVW8dQgw1cs1LLGbSJUJtA/8CvpPEWlb1a+uItI48KvthyLysYgMqOpkInJ34JhPReThwHMDReT1wIrAhMCvum2AR4DLAte+rLpriMjVIvJK4C9koYj8I/B8AxGZKCKfi8hnInJH4PmJIjIo8HPPwLk+C5y7UeD51SLysIgsCbxW1V9wE3uygVtU9TTgT7jOl/uJyPG4LphveDlZxK88NXFrOHBy4FdQRKRtNcf9BXhDVYeJyNHAYhF5PdBvhcB7+wAdgK64VY4zReRcVX1VRC7BtQTuCzyoqmtE5K9AqqreHHj/36u6RuD0XYBTgTLgKxEZDTQHWqnqyYH3H10xYBE5FPfreE9VXSEik4AbcA3PwPWL+U2gFPUn4LoD+i9nooqIHAF0B6bJz43uGlU6bDAwXVX3ejmnjdhNtOsDDBeRT4AFwKFAmyqO6QN8DCzBNZvqEHjtFlxv8DJVre7X3JquMV9Vv1fVncBy4Hhcz5MTRGS0iPQFKncsTAZWVej1kwOcW+H1fXXUj3B9vk1sOwj4bl8tPfDVsdIxg/FYhgEbsZvoJ8AlqvpVLceMUNXnq3itFa4vSwsROaianuBVXkNEzsCN1PfZCzRU1a0icgquNe1NQDowrNL5arLvnHuxv6MxT1V/EJFVInKpqk4TN2xPUdWlACKSjNtU532v57QRu4lU24AjPRyXD9wS+MuAiJxazTHDAr/yIiKtRKS5iDQEJgC/x3VRvLOaa3u5xn4i0hS3G9DLwAPAbyod8iXQVkSSAo+vBN7ycK8mBgSahr0PJItIkYhcC1wOXCsiS4FlQMXPioYAuXoAjb1sNGAikqpuFpF3Ax+YzgGyqjn0b7ja9KeBxLsauLjSuf4rIh1xG1cAbAeuAP4ILFTVhYEyy4ciMgt4k59LLyO8XKOSVsAEEdk3cLq3Ujw7ReQaXE21Ia7Nbcx0sjQ1U9Uh1bxU5RRIVX3oQK9h3R2NMSbGWCnGGGNijCV2Y4yJMZbYjTEmxlhiN8aYGGOJ3RhjYowldmOMiTGW2I0xJsb8f7awPGfHcCKwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tile_area = 256*256 \n",
    "tile_area\n",
    "z_level = np.linspace(0,20,21)\n",
    "tiles = 2**(2*(z_level))\n",
    "km_pixel = ((156412)/(np.sqrt(tiles)))/1000\n",
    "km2_pixel = ((156412)/(tiles))/1000\n",
    "exte = np.round(tile_area * km2_pixel)\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.ylabel(\"km/pixel (min grid area)\")\n",
    "plt.xlabel(\"tile extension\")\n",
    "\n",
    "plt.plot(  exte, km_pixel, 'r')\n",
    "\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "id": "034d08b6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD9CAYAAABHnDf0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnvUlEQVR4nO3de3xU9Z3/8ddHQEAU1BIwGARsYrgoRh3RVlupGMVLk223i7ja0qplW211t60We1nbbrektRetsmtjUeNqTbEXY9WCiKL9WZUOeAWF1FKFECEI3rhZwuf3x5nBGGeSSTIzZy7v5+NxHjNzzpk5n5wMn3z5nu/5fM3dERGRwrJP2AGIiEj6KbmLiBQgJXcRkQKk5C4iUoCU3EVEClD/sAOI6fWQnenTp7Nw4cJ0xiIiki8s2Ya8b7lv3rw57BBERHJO3id3ERF5v4wmdzO72cw2mdnzndZ/2cxWm9lKM/tRfP3cuXMpLy+nsrKSRYsWZTI0EZGCluk+91uBG4Db4ivM7GNALTDZ3XeZ2QjgilWrVtHY2MjKlSvZsGEDp512GmvWrKFfv34ZDlFEpPBktOXu7o8CWzqt/iJQ5+67YvtsAmhqamLmzJkMHDiQcePGUV5ezrJlyzIZnohIwQqjz/0I4CNm9qSZPWJmxwO0tLQwevTovTuVlZXR0tKS8APq6+uJRCJEIhHa2tqyErSISD4JI7n3Bw4CTgSuABa4O4kKmJklHuUze/ZsotEo0WiUkpKSTMYqIpKXwkju64HfeWAZsGfz5s2UlZWxbt26d3dav55Ro0aFEJ6ISP4LI7nfDZwKYGZHAPsOHz6cmpoaGhsb2bVrF2vXrqW5uZkpU6Z0/UkPPABJum5ERIpZRkfLmNmdwFRguJmtB64GbgZujg2PfAeYZWZLJk2axIwZM5g4cSL9+/dn3rx53Y+UefJJePVV2L4d9tsvkz+KiEhesRyZrKN3QSxYQOTcc4k+9RRUVaU3IhGR3Feg5QfGjw8eV68ONw4RkRyT38m9oiJ4fPHFcOMQEckx+Z3cBw+GffdVy11EpJP8Tu4Agwap5S4i0klhJPc1ayA3LgyLiOSEwkju27ZpvLuISAeFkdxBXTMiIh0UTnLXRVURkb3yP7kPGAD776+Wu4hIB/mf3CG4mUktdxGRvQojuVdWKrmLiHRQOMn9lVeCUTMiIlIgyT1eY6a5Odw4RERyRGEk98rK4FEXVUVEgEJJ7hUVYKZ+dxGRmMJI7oMHw5gxarmLiMRkNLmb2c1mtik261LnbV8zMzez4fF1c+fOpby8nMrKShYtWtSzg2nEjIjIXpluud8KTO+80sxGA9XAK/F1q1atorGxkZUrV7Jw4UIuueQS2tvbUz9SfKz7nj19j1pEJM9lNLm7+6PAlgSbfgZcSYfp9Zqampg5cyYDBw5k3LhxlJeXs2zZstQPVlkZzKWqAmIiItnvczezGqDF3Z/puL6lpYXRo0fvfV1WVkZLkkRdX19PJBIhEonQ1tYWrNSUeyIie2U1uZvZfsA3gf/svC3RRN1mied+nT17NtFolGg0SklJSbBSwyFFRPbqn+XjfRAYBzwTS9xlwIpXX32VsrIy1q1bt3fH9evXM2rUqNQ/ubQUDjhALXcREbLccnf359x9hLuPdfexwHrg2EMOOYSamhoaGxvZtWsXa9eupbm5mSlTpqT+4WYaMSMiEpPRlruZ3QlMBYab2Xrganefn2jfSZMmMWPGDCZOnEj//v2ZN28e/fr169kBKyvh0Uf7GraISN6zRH3dIeh1EJFIhGg0Grz4/vfh29+Gt9+GIUPSFZuISK5KfGGSQrlDNS5+UXXNmnDjEBEJWWEldw2HFBEBCi25l5cHF1Y1HFJEilxhJfd4ATG13EWkyBVWcgfNpyoiQiEm9/hYdxUQE5EiVnjJffx4FRATkaJXeMldNWZERAo4uavfXUSKWOEldxUQExEpwOQeLyCmbhkRKWKFl9xBwyFFpOgVZnKvrIR162DbtrAjEREJRWEm93iNGRUQE5EilVI9dzOLAB8BRgE7gOeBB9090eTX4es4HPKYY8KNRUQkBF223M3ss2a2ArgKGAysBjYBJwOLzazBzA7LfJg9FC8gpn53ESlS3bXchwAnufuORBvNrAqoAF5Jsv1m4Bxgk7sfGVt3DfBx4B3gJeBz8QlD5s6dy/z58+nXrx8///nPOeOMM3r+E0FQQGzsWCV3ESlaXbbc3X1essQe2/60uy/p4iNuBaZ3WrcYONLdJwNrCP5XwKpVq2hsbGTlypUsXLiQSy65hPb29pR+iIQ0HFJEilhKF1TN7AgzW2Jmz8deTzazb3X3Pnd/FNjSad0D7r479vIJoAygqamJmTNnMnDgQMaNG0d5eTnLli3r0Q/zHuPHBxdUVUBMRIpQqqNlbiJoYf8DwN2fBWam4fgXAn8EaGlpYfTo0Xs3lJWV0ZKk+Fd9fT2RSIRIJEJbW1viT66sDAqIrV+fhjBFRPJLqsl9P3fv3IzenXDPFJnZN2OfcQdAoom6zRLP/Tp79myi0SjRaJSSkpLEB1CNGREpYqkm981m9kHAAczsU0Brbw9qZrMILrSe77GsXlZWxrp16/bus379ekaNGtXbQ7w71l397iJShFJN7pcCvwDGm1kL8O/AF3tzQDObDnwdqHH37fH1NTU1NDY2smvXLtauXUtzczNTpkzpzSEChxyiAmIiUrRSuonJ3f8GnGZmQ4B93P2tVN5nZncCU4HhZrYeuJqg734gwTh5gCfcnUmTJjFjxgwmTpxI//79mTdvHv369evFj7T34KoxIyJFyxL1db9vJ7N24Brgqng3ipmtcPdj0xRH90EkEYlEiEajiTd++tOwdGlQZ0ZEpPAkvjBJ6t0yK2P7PmBmB3f3oTlj/PhgtMzbb4cdiYhIVqWa3He7+5UEQyL/ZGbH0YfWdtbER8yogJiIFJlUk7sBuPsCYAZwC3B4poJKGw2HFJEildIFVeDi+BN3X2lmJwP/lJGI0qmiQgXERKQodZnczexUd38IGGNmYzptzv2O7EGDggJiGusuIkWmu5b7KcBDBFUcO3Pgd2mPKN00HFJEilCXyd3dr449fi474WRAZWUwHHLPHtinMCeeEhHpLNWqkJeb2VAL/NLMVpjZ6ZkOLi0qK2HHDhUQE5GikmpT9kJ3fxM4HRgBfA6oy1hU6aQaMyJShHo0FBI4C7jF3Z8hH25iAg2HFJGilGpyX25mDxAk90VmdgCQH7NgHHIIDB2q5C4iRSXVce4XAVXA39x9u5l9gKBrBgAzm+TuKzMQX9+Zaco9ESk6KbXc3X2Pu69w99djr1+LzcYU93+ZCC5tNBxSRIpMusYG5nb/e2WlCoiJSFFJV3LP7SJiKiAmIkWmOO7qiQ+HVNeMiBSJdCX3dxKtNLObzWyTmT3fYd3BZrbYzJpjjwfFt82dO5fy8nIqKytZtGhRmkIDysuDC6u6qCoiRaK7wmFdzrTk7itijycm2eVW4Abgtg7r5gBL3L3OzObEXrNq1SoaGxtZuXIlGzZs4LTTTmPNmjV9m2ovbtAgGDdOLXcRKRrdDYX8SexxEBAB4jcvTQaeBE7u6s3u/qiZje20upZgXlWABmApQFNTEzNnzmTgwIGMGzeO8vJyli1bxoc+9KEUf5RuaDikiBSRLrtl3P1j7v4x4GXgWHePuPtxwDHAX3t5zJHu3hr7/FaCcga0tLQwevTovTuVlZXR0tLSy0MkUFkZXFDdkx/3XomI9EWqfe7j3f25+At3f57gpqa0STRRt1niEZb19fVEIhEikQhtbW2pHWD8+KCAmCbLFpEikGpyfyFWDXKqmZ1iZjcBL/TymBvNrBQg9rgJgpb6ug6Jd/369YwaNSrhB8yePZtoNEo0GqWkpCS1o6rGjIgUkVST++eAlcDlwL8Dq+hQfqCH7gFmxZ7PApoAampqaGxsZNeuXaxdu5bm5mamTJnSy0MkoOGQIlJEUqot4+47gZ/FlpSZ2Z0EF0+Hm9l64GqCUsELzOwi4BXgX4ArJk2axIwZM5g4cSL9+/dn3rx56RkpEzdyZFBATBdVRaQIWKK+7r0bzRa4+wwze44Ed6G6++Q0xdHrO1wjkQjRaDS1nU84AQ44AB58sLeHExHJJUlLv3TXcr889nhO+mIJUWUlPPRQ2FGIiGRcd0MhW82sHzDf3V/uvGQpxvSprISWFhUQE5GC1+0FVXdvB7ab2bAsxJNZ8YuqKiAmIgUu1ck6dgLPmdliYFt8pbtflpGoMiU+HPLFF+HYLisriIjktVST+32xJb+Vl8M++2g4pIgUvFSHQjZkOpCsGDQIxo7VcEgRKXgpJfckQyHfAKLA9939tXQHljGVlWq5i0jBS7Vb5o9AO/Cr2OuZBOMr3yAo6/vxtEeWKePHw9KlQQGxfYpjrhIRKT6pJveT3P2kDq+fM7PH3P0kM7sgE4FlTGXluwXExowJOxoRkYxItem6v5mdEH9hZlOA/WMvd6c9qkxSjRkRKQKpttwvBm42s3hCfwu4yMyGAHMzElmmdBwOefrp4cYiIpIhqY6W+QtwVOxGJnP31ztsXmBms/JmRM3IkTBsmFruIlLQenRF0d3f6JTY4y5PsC43mWnKPREpeOkaLpK0MllO0nBIESlw6UruvS7ZG4rx44MCYm+9FXYkIiIZUbwtd1ABMREpWOlK7o+l6XOyQ8MhRaTApZTczWygmf2rmX3DzP4zvsS3u/uXenpgM/sPM1tpZs+fd9557Ny5ky1btlBdXU1FRQXV1dVs3bq1px+bmngBMV1UFZEClWrLvQmoJbhhaVuHpVfM7FDgMiDi7ke2t7fT2NhIXV0d06ZNo7m5mWnTplFXV9fbQ3Rt4MCggJha7iJSoFJN7mXufq67/8jdfxJf+njs/sBgM+u/fft2Ro0aRVNTE7NmzQJg1qxZ3H333X08RBfGj1fLXUQKVqrJ/c9mdlS6DuruLcCPgVeA1mHDhnH66aezceNGSktLASgtLWXTpk0J319fX08kEiESidDW1ta7ICorobk5KCAmIlJgUk3uJwPLzWy1mT1rZs+Z2bO9PaiZHUTQzTMOGLVt2zZuv/32lN8/e/ZsotEo0WiUkpKS3gUxfnxQQOyVV3r3fhGRHJZqbZkz03zc04C17t4GcNttt/HnP/+ZkSNH0traSmlpKa2trYwYMSLNh+0gEgkeH3kk6H8XESkgXbbczWxo7OlbSZbeegU40cz2MzNbsmQJEyZMoKamhoaGoERNQ0MDtbW1fThEN445Bg49FJqaMncMEZGQdNdy/xVwDrCc4C7UjjcrOXB4bw7q7k+a2W+AFcDuPXv2MHv2bN5++21mzJjB/PnzOeyww7jrrrt68/GpMYOaGmhoCLpnBg/O3LFERLLM3HOickCvg4hEIkSj0d69edEimD4d7r0Xzj67tyGIiIQlaXWAVPvcMbPJwNiO73H33/UprLBNnQoHHBB0zSi5i0gBSXWC7JuBycBKID520IH8Tu4DB8KZZ8If/qA5VUWkoKTacj/R3SdmNJKw1NbCggWwbBmceGLY0YiIpEWqTdXHzawwk/uZZ0K/fho1IyIFJdXk3kCQ4NNyE1NOOeggOOUUJXcRKSipdsvcDHwaeI53+9wLR20tXH55UI6goiLsaERE+izVlvsr7n6Pu69195fjS0Yjy6b4zVL33BNuHCIiaZJqcn/RzH5lZueZ2SfjS0Yjy6YxY+Doo9U1IyIFI9XkPhjYBZwOfDy2nJOpoEJRUwOPPQabN4cdiYhIn6XU5+7un8t0IKGrrYX/+q/gbtXPfjbsaERE+qS7wmHfMrODu9h+qpkVRgv+2GOhrExdMyJSELpruT8H/MHMdhIU+WoDBgEVQBXwIPCDTAaYNfFCYrfeqkJiIpL3umy5u3uTu58EfIGg9EA/4E3gdmCKu/9HvCZ7Qaithe3bYcmSsCMREemTVPvcm4HmDMcSvo6FxM4pjN4mESlOqpTV0b77vreQmIhInlJy76y2FjZuhCefDDsSEZFeCy25m9mBZvYbM3txwoQJPP7442zZsoXq6moqKiqorq5m69at2Q/srLOgf3/drSoiea3LmZjM7Hq6mCXJ3S/r9YHNGoA/ufsv33nnHd++fTs/+MEPOPjgg5kzZw51dXVs3bqVH/7wh11+Tp9mYkrmtNNgwwZYtSq9nysikl5JZ2LqLrnP6upT3b2hV9EEE28/AxzuQQAOUFlZydKlSyktLaW1tZWpU6eyevXqLj8rI8n9+uvhsstgzRoVEhORXNa75P6+nc2GuPu2PkdjVgXUA6uAoy+66KKq6667jkMPPZTXX399734HHXRQwq6Z+vp66uvrAWhra+Pll9Ncw+zll2HsWLjmGvja19L72SIi6ZM0uafU525mHzKzVcALsddHm9n/9CGg/sCxwP+6+zFDhgyhrq4u5TfPnj2baDRKNBqlpKSkD2EkoUJiIpLnUr2gei1wBvAagLs/A3y0D8ddD6x39ycBPvWpT7FixQpGjhxJa2srAK2trYwYMaIPh+ij2lr485+hrXDu0RKR4pHyaBl3X9dpVXtvD+rurwLrzKwSYMmSJUycOJGamhoaGoJu/IaGBmrjddbDUFsbjHW/777wYhAR6aVUZ2JaZ2YfBtzM9gUuI9ZF0wdfBu4ws31ra2u55ZZb2LNnDzNmzGD+/Pkcdthh3HXXXX08RB8ccwyMHh10zahKpIjkmVST+xeA64BDCbpUHgAu7cuB3f1pIBJ/GV+/JFfqusQLid1yiwqJiUjeSbVbZh93P9/dR7r7CHe/ABieycByQk1NUEjswQfDjkREpEdSTe5/MrMZ8Rdm9lXg95kJKYdMnQpDh2rUjIjknVS7ZaYC9Wb2L8BIgv72KZkKKmd0LCTW3g79+oUdkYhISlJqubt7K7AQ+BAwFrjN3d/OYFy5o7YWNm2CZcvCjkREJGWp3sS0GDgBOBI4C/iZmf04k4HljDPPDAqJqWtGRPJIqn3u89z9M+7+urs/D3wYeCODceWOAw+EU05RcheRvJJqt8zdnV7vdvf/ykhEuai2Fl58MSgkJiKSB7pM7mb2/2KPb5nZm7HH+PJmdkLMATU1waNa7yKSJ7qbIPvk2OMB7j409hhfhmYnxBwwZgxUVWkCDxHJGynXljGzY83sMjP7spkdk8mgcpIKiYlIHkl1tMx/Ag3ABwjuTL3VzL6VycByTryQ2L33hh2JiEi3Upqsw8xeAI5x952x14OBFe4+IU1xpD5jSCcZmYkpEfege+bYY+HuuzN/PBGR7vVtsg7g78CgDq8HAi/1IaD8Ey8k9sADQb0ZEZEclmpy3wWsNLNbzewW4HngbTP7uZn9PHPh5Zja2qBCpAqJiUiOS7W2zO95b6GwpekPJQ+cckpQSOyee94dHikikoO6TO5mVg/8Efidu7+VnZBy2L77wllnqZCYiOS87rplbgaOBu43syVm9nUzOzpdBzezfmb21DnnnAPAli1bqK6upqKigurqarZu3ZquQ6VPTU1QSOzJJ8OOREQkqe5uYnrC3b/j7h8BZgCvAF81s6fN7OaONd576XI6TNdXV1fHtGnTaG5uZtq0adTV1fXx4zNAhcREJA/0ZILs19z9zlgBsSpgHlDR2wObWRlwNvDL+LqmpiZmzZoFwKxZs7g7F4ccHnhgMImHkruI5LCULqia2YHAZwhque99j7tf1odjXwtcCRwQX7Fx40ZKS0sBKC0tZdOmTQnfWF9fT319PQBtYdwxWlsLX/4yrF4NlZXZP76ISDdSbbnfT5DYnwOWd1h6xczOATa5e68+Y/bs2USjUaLRKCUlJb0No/dqa4Nx7//zP9k/tohIClIdCjnI3b+SxuOeBNSY2VnAoIceeogLLriAkSNH0traSmlpKa2trYwYMSKNh0yj0aPhC1+AG26Az34Wjim+UjsikttSbbn/n5l93sxKzezg+NLbg7r7Ve5e5u5jgZmnnnoqt99+OzU1NTQ0NADQ0NBAbW1tbw+ReT/4AQwfDv/2b8GwSBGRHJJqcn8HuAZ4nHe7ZNJe0GXOnDksXryYiooKFi9ezJw5c9J9iPQ58ED46U/hL3+BWP+/iEiuSLVw2EvACe6+OUNx5H7hsETc4bTTYPny4OLqyJHhxCEixarPhcNWAqqW1Vn8ouqOHfDVr4YdjYjIXqkm93bgaTP7RbxYWFEVDOtKZSV8/etwxx2wZEnY0YiIAKl3y8xKtN7dG9IUR352y8Tt2AFHHRXUmnn2WRg4MNx4RKRY9Llb5nl3b+i4AFvSE1sBGDwY5s2DNWvgRz8KOxoRkZST+01mdlT8hZmdBxTXNHvdOeMMmDED/vu/4a9/DTsaESlyqSb3TwENZjbBzD4PXAKcnrmw8tTPfhaUBb700mAkjYhISFJK7u7+N2Am8FuCRH+6u7+RycDy0qhRQcv9gQfgrrvCjkZEiliXF1TN7Dnee7FzBPAGwbR7uPvkNMWR3xdUO2pvhylToLUVXnwxmLlJRCQzkl5Q7a62zDlpDqTw9esHN94IJ5wA3/42XHdd2BGJSBHqMrm7+8vZCqSgHH88XHJJUFjsM5+B444LOyIRKTIpT9YhPfT970NJSVA9UoXFRCTLlNwz5cADg9Ez0Sj84hdhRyMiRUbJPZNmzgwKi111Fbz6atjRiEgRUXLPJLPgztWdO+Er6ZzrRESka0rumXbEEUHL/c474cEHw45GRIpESoXDsqBwxrknsnNnUFjMLCgsNmhQ2BGJSGHoc+GwtDKz0Wb2sJm9YGYrr4uNBd+yZQvV1dVUVFRQXV3N1q1bwwgv/QYNCuq+NzfDD38YdjQiUgTC6pbZDXzV3ScAJ86bN49Vq1ZRV1fHtGnTaG5uZtq0adTV1YUUXgZUVwcXWOfODZK8iEgGhZLc3b3V3VfEnr81YcIEWlpaaGpqYtasoHT8rFmzuPvuu8MIL3N++tOg1rsKi4lIhoV+QdXMxj711FOccMIJbNy4kdLSUgBKS0vZtGlTyNGlWWlpUFhs8WJoSNc8JyIi7xdqcjez/YHfXnvttQztQYGt+vp6IpEIkUiEtra2zAWYCV/8Ipx8Mlx4YdBFoxa8iGRAaKNlzGwAcC+wyN1/AlBZWcnSpUspLS2ltbWVqVOnsnr16i4/Jy9Gy3S2fTtcfHEwPPLcc2H+fBgyJOyoRCT/5NxoGQPmAy+4+0/j62tqamiIdVc0NDRQW1sbRniZt99+wYTaP/oRLFgQtORfVo02EUmfsLplTgI+DZxqZk9XVVVx//33M2fOHBYvXkxFRQWLFy9mzpw5IYWXBWZwxRVw332wdi1EIvDII2FHJSIFQjcx5YI1a6C2Nph79ec/DypJWtL/bYmIxOVWt4x0csQR8MQTwSTbl1wSJPd33gk7KhHJY0ruuWLYMGhqgm98A+rr4dRTYePGsKMSkTyl5J5L+vULxsH/+tfw1FNBP3y+dzmJSCiU3HPRjBnw2GNBsv/IR4KRNSIiPaDknquqquAvfwkm2r7ggmBkjabrE5EUKbnnspKSoFTBpZfCj38MZ58NhVIpU0QySsk91w0YADfcADfdBA89BFOmwPPPhx2ViOQ4Jfd8cfHF8PDD8NZbMHkynHMO3HuvumpEJCEl93xy0knw9NPwrW/BihXw8Y/DuHHw/e9Da2vY0YlIDlFyzzeHHALf+15Qi+a3v4Xx4+Hb34bRo+Gf/znoo9+zJ+woRSRkSu75asAA+OQn4YEHgpmdvvKVoDbN6acHd7xecw3kWzlkEUkbJfdCUF4eVJhsaQnGxB96KFx5JZSVwfnnw5/+pLrxIkVGyb2QDBwI//qvQQt+5cqgRs1998FHPwpHHgnXXw+vvx52lCKSBaoKWei2bw/KGdx4IyxbBv37B4n+uOOC8gaRCBx1VPCHQUTyTdKqkEruxeSpp+A3vwnq1USjsGVLsH7AgGB4ZTzZRyIwaVKwXkRymZK7dOIejLiJJ/r48sYbwfaBA4MSCB0T/vjxQctfRHKFkrukYM8e+Nvf3pvsly+Ht98Otg8aFAy5HDUqWEpL333ecdF8sCLZkj/JfeHChVx++eW0t7dz8cUXdzvVnpJ7hu3ZEwy1jEaDbp2WFtiwIVhaWmDHjve/Z+jQxEn/oIPggAMSL/vvr/8ViPRcfiT39vZ2jjjiCBYvXkxZWRnHH388d955JxMnTkz6RiX3ELnDm2++m+y7WlKZWWrw4OSJf7/9gq6iffcNHjs+7+5xwICgfHK/fsEfkPjzRK8Trdtnn2Axe/9zM02JKGFK+uXLqabSsmXLKC8v5/DDDwdg5syZNDU1dZncJURmwQxSw4bBhAnJ93MPLt6+/npQG6cny8aNwdyyO3bArl3BH4n4Yy7V1UmU9Dsn//jzZOuS7RNf35PHzs/TsS5VYb23r8I69s9+FtSKSrOcaLlPnz7dN2/ezNatW3nzzTcZM2YMAK+99hrbtm3jsMMOe8/+bW1tbN68GYBdu3ZRVVWV7ZC71dbWRklJSdhhvE/BxbVnT/DHwz35844L9Oj5tm3bGDJkyHtvAuvp80Sve7s+ZufOnQwaNCilfXusj5+3a9cuBubg0NpdO3cyMH7OcsjrAwZwYFlZr967fPnyRe4+PeFGd8+Fxd3dFyxY4BdddFH8pd92223+pS99ybuy3377dbk9LMcdd1zYISSkuHpGcfVcrsZWoHElzas5dYdqWVkZ69at2/t6/fr1jBo1KsSIRETyU04l9+OPP57m5mbWrl3LO++8Q2NjIzU1NWGHJSKSd3Lqgmr//v254YYbOOOMM2hvb+fCCy9k0qRJXb5n+PDhWYquZ2bPnh12CAkprp5RXD2Xq7EVW1w5cUEV3cQkItIbSYf45FS3jIiIpIeSu4hIAcqb5L5w4UIqKyspLy+nrq7ufdvdncsuu4zy8nImT57MihUrMh7TunXr+NjHPsaECROYNGkS11133fv2Wbp0KcOGDaOqqoqqqiq+973vZTwugLFjx3LUUUdRVVVFJBJ53/Ywztfq1av3noeqqiqGDh3Ktdde+559snW+LrzwQkaMGMGRRx65d92WLVuorq6moqKC6upqtm7dmvC93X0X0x3XFVdcwfjx45k8eTKf+MQneD1JTf7ufufpjus73/kOhx566N7f1f3335/wvdk+X+eee+7emMaOHZv0PphMnq9kuSGr37GuxklmcenS7t27/fDDD/eXXnrJd+3a5ZMnT/aVK1e6+7tjRO+77z6fPn2679mzxx9//HGfMmVKb8aM9siGDRt8+fLl7u7+5ptvekVFxd644h5++GE/++yzMx5LZ2PGjPG2trak28M4Xx3t3r3bR44c6X//+9/fsz5b5+uRRx7x5cuX+6RJk/auu+KKK3zu3Lnu7j537ly/8sor3/e+rr6LmYpr0aJF/o9//MPd3a+88sqEcbl3/ztPd1xXX321X3PNNV2+L4zz1dFXvvIV/+53v5twWybPV7LckIHvWNK8GnZST2kBPgQs6vD6KuCq2POFscdfAOd12Gc1UJrlOJuA6k7rpgL3hnDO/g4M72J7qOcLOB14LMH6rJ0vYCzwfKJzAJQCqxO8J+l3MVNxddr2CeCO3vzOM3C+vgN8rZv3hHa+CC42rgMqwjhfnY7VBFRn8zuWL90yhxL8kuLWx9bh7956m3SfbDCzscAxwJMJNn/IzJ4xsz+aWddjO9PHgQfMbLmZJRprFer5AmYCdybZFsb5Ahjp7q0AsccRCfYJ+7xdCPwxybbufueZ8CUze9bMbjazgxJsD/N8fQTY6O7NSbZn5Xx1yg1Z+47lS3JPNNyn8/DJVPbJCDPbH/gt8O/u/manzSuAMe5+NHA9cHc2YgJOcvdjgTOBS83so522h3m+9gVqgLsSbA7rfKUqzPP2TWA3cEeSXbr7nafb/wIfBKqAVuAnCfYJ7XwB55G8AQFZOF/d5Iakb0uwrsfnLF+S+3pgdIfXZcCGXuyTdmY2gOCXd4e7/67zdnd/093fjj2/HxhgZhm/88rdN8QeNwG/B6Z02iWU8xVzJrDC3Td23hDW+YrZaGalALHHTQn2Cet7Ngs4BzjfY/9X7yyF33lauftGd2939z3ATUmOF9b56g98Evh1sn0yfb6S5IasfcfyJbn/Bagws3GxVt9M4J5O+9wDfMYCJwJvxP/7kylmZsB84AV3/2mSfQ6J7YeZTSE4569lOK4hZnZA/DlB//bznXbL+vnqIGmLKozz1cE9wKzY81kE/aSdpfJdTCszmw58Hahx9+1J9knld57uuEo7vPxEkuNl/XzFnAa86O7rE23M9PnqIjdk7zuWjYsJabogcRawBngJ+GZs3ReAL/i7F0/mxbY/B0SyENPJBP9dehZ4Orac1SmuLwErgWeAJ4APZyGuw2PHeyZ27Jw4X7Hj7keQrId1WJf180Xwx6UV+AdBS+ki4APAEqA59nhwbN9RwP1dfRczHNdfCfpg49+xGzvHlex3nuG4/i/23XmWIPmUdo4rjPMVW39r/DvVYd9snq9kuSFr37FcKT8gIiJplC/dMiIi0gNK7iIiBUjJXUSkACm5i4gUICV3EZECpOQuIlKAlNwl55jZ39N9V2qyz4zdxPWQmQ1NwzHazezpDsucbvb/Rh+P909mNrHD6x+b2al9+UwpHEruUuzOAp7xWN2PJMWvUrXD3as6LN0V4u5Tcgf+CZjY4fX1QJd/UKR4KLlLVpjZFzq0aNea2cMpvu8CM1sWe98vzKyfmX3RzH7UYZ/Pmtn1yfbv5hDn895bwKNm9iszOzVeBqEvzGyYma02s8rY6zvN7PNmVgcMjsV5R1exm9nbZvbfFlTKfMLMRprZhwmKr10T2/+D7v4y8AEzO6SvcUv+U3KXrHD3G929Cjie4DbxhLV4OjKzCcC5BNX7qoB2gmT8G4KiUHHnAr/uYv+unAQs7/D6COBXBGUQVpnZN8xsVHexxsSTdXw5193fiH3WrWY2EzjI3W9y9zm829I/v5vYhwBPeFAp81Hg8+7+Z4Jb/q+IfcZLsX1XxH4mKXL9ww5Ais51wEPu/ocU9p0GHAf8JdaIHgxscvc2M/tbrOBZM1AJPAZcmmj/bo5xsLu/FX/h7u3AvcC9ZlYCzAVeMbMPu/uybj5rRywxv4e7LzazfyGo5XN0T37W2LZ3YjFB8IeouosYNhHUKZEip+QuWWNmnwXGELRkU3oL0ODuVyXY9mtgBvAi8Ht391g3SrL9k9ltZvt4ULY2Hucwglb05wgKUl1EUACq488yGoj/gbrR3W9M+kOY7QNMAHYABxP8z+V9u3UR+z/83SJQ7XT973ZQ7DhS5NQtI1lhZscBXwMu6JhIu7EE+JSZjYh9xsFmNia27XcEFxTP492a3V3tn8xqggqB8ThvJ+jaOBz4jLt/1N0b3H1nxze5+7oOF06TJvaY/wBeiMV6swV1vgH+0eF5b2J/Czig07ojyHCpX8kPSu6SLV8iaLU+HOuP/mV3b3D3VcC3CKZCexZYTDDvJO6+FVhFMGvTsu7278J9BPO2xi0AKt19jiefni2Zzn3udWZ2BHAx8FV3/xNBn/m3YvvXA8+a2R29jL0RuMLMnjKzD8b+UJQD0R7GLQVIJX+lqFkw4cRt7t5VP3ZeMLNPAMe6+7fDjkXCp5a7FDUPZp+6KR03MeWA/iSex1SKkFruEhozexIY2Gn1p939uTDiESkkSu4iIgVI3TIiIgVIyV1EpAApuYuIFCAldxGRAvT/AQIo8Q9skQD8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.xlabel(\"z_level (> -Extent)\")\n",
    "plt.ylabel(\"km/pixel (min_grid_size)\")\n",
    "\n",
    "# plot the function\n",
    "plt.plot(z_level, km_pixel, 'r')\n",
    "\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "id": "190a9432",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD8CAYAAABuHP8oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfU0lEQVR4nO3deXxU9b3/8deHRaRuLALGBVGJSlFwiaUovUZyU6Ei3FJLtfordcOt1dZLFe6tlbZqgxsKuBBFiYpoQSpUBUQqahUvq1gBlSouYEqABEGKmITv74/vgAhZJmTOOXNm3s/HI4/JnDkz53MekHdOvue7mHMOERGJnyZRFyAiIntHAS4iElMKcBGRmFKAi4jElAJcRCSmmoV8vL3u8tKnTx9mzpyZylpEROLAanshNlfg69evj7oEEZG0EpsAFxGRb1KAi4jElAJcRCSmFOAiIjGlABcRiSkFuIhITCnARURiKuyBPCIimePZZ2HatPr3a9YMHnoo5YdXgIuI7I3ly2HQIDjgANh//7r3bdEikBIU4CIiDbV9O1x+uQ/vFSugfftIygilDdzMmprZkn79+gFQXl5OYWEhubm5FBYWUlFREUYZIiKpMW4cvPEG3HVXZOEN4d3EvA5YseNJUVERBQUFrFy5koKCAoqKikIqQ0SkkdasgRtvhIICGDw40lICb0Ixs8OBc4BbgQsApk2bxty5cwEYPHgw+fn5jBw5MuhSRERqt2UL/Pzn8OGHde9XVgZVVf4q3GqdKDAUYbSB3wPcABywY8PatWvJyckBICcnh7KyshrfWFxcTHFxMQDr1q0Luk4RyWYjRsCUKdC3LzRtWvt+hx4KF18MxxwTWmm1CTTAzawfUOacW2Rm+Q19/5AhQxgyZAgAeXl5qS1ORGSHRYvg7rthyBB/ZR0TQbeBnwH0N7OPgKf+9re/cdFFF9GhQwdKS0sBKC0tpX2ENwFEJMtVVfkeJR06QMyacgMNcOfccOfc4c65TsD5vXv35oknnqB///6UlJQAUFJSwoABA4IsQ0SkdqNGwZIlMGYMtGoVdTUNEkk/8GHDhjFo0CDGjx9Px44dmTx5chRliEimmTTJX01v25b8e6qqYMAAGDgwuLoCYs7t9TKVe2OvD5aXl8fChQtTWYuIZJLSUujSxd9c7NMn+fe1bAlXXQVt2wZXW+PU2tVFIzFFJDNcey18+SU89RTk5kZdTSgU4CISf9On+y6At96aNeENCnARSWeVlfDPf9a/z9VXw4knwm9+E05daUIBLiLpqboazjwT5s2rf18zeOYZaN48+LrSiAJcRNLT6NE+vH//ezjuuLr3zc2FU04Jp640ogAXkfTz0Ufw29/COefATTdFPudIutKSaiKSXpyDK6+EJk3g/vsV3nXQFbiIBOvNN+GTT5Lff/lymDUL7r0XOnYMrq4MoAAXkeA8/zwkFnJpkF694JprUl9PhlGAi0gwNm/2Ixy//W34858b1hSSm1v3lK4CKMBFJCi//S2sXg2vvw5du0ZdTUbSTUwRSb3/+z8/u9/VV0PPnlFXk7F0BS4iyVu/Hh580M85UpcpU/zKNbfdFk5dWUoBLiLJcc6vGfn88/W3T++/v59U6sADQyktWynARSQ5f/6zD++77oLrr4+6GkFt4CKSjPJyP11rXp5/lLSgK3ARqd/QobBhA7z4IjRTbKQL/UuIZLMpU+Cxx+rep6oKZsyAYcOge/dw6pKkKMBFstX778NFF0G7dv6rLj/+Mfzud+HUJUlTgItko+3bYcgQ2HdfmD8fcnKirkj2ggJcJBs98gi88goUFyu8Y0y9UESyzb/+5ZceO/NMuPTSqKuRRtAVuEim2LoVzjsPli2re7/Nm/2+xcV+zm2JLQW4SKb44x/hhRfg/POhRYu69x04EI49Npy6JDAKcJFMsHQp3H47XHyxb9+WrKC/n0TirroaLr8c2raFO++MuhoJka7ARdLZl19CZWXd+xQXw4IFMGkStGkTTl2SFhTgIulqxgzfVl3f1K0AP/gB/OQnwdckaUUBLpKONm3yA206dYLLLqt733328SMqtXp71lGAi6Sj//1fWLMG3ngDvvvdqKuRNKWbmCLpZt48uO8++MUvFN5SJ12Bi4SlshIWLvTzkNTGObjySjj8cLj11vBqk1hSgIuEwTl/Q/K555Lb/69/hQMOCLYmiT0FuEgYnnrKh/ewYdC7d937HnIInHhiOHVJrCnARYK2YQNcdx185ztwyy31LwgskiQFuEjQ/vu/oaICXnpJ4S0ppQAXSYV33oHXXttz+7p1UFICw4dDt27h1yUZTQEu0lhr1sDpp/tpWmvSvTvcdFO4NUlWUICLNIZzcM01fuHfxYvh0EP33KdNG2jePPzaJOMFOpDHzPY1s/lmttTMlt18880AlJeXU1hYSG5uLoWFhVRUVARZhkhwpk6FadNgxAg4+WTo0GHPL4W3BCTokZjbgN7Oue7ASTNnzuTNN9+kqKiIgoICVq5cSUFBAUVFRQGXIRKAjRv9aMmTT4brr4+6GslCgTahOOcc8EXiafPKykrMjGnTpjF37lwABg8eTH5+PiNHjgyyFJHkVVX5xRE+/bTu/ZYtg7Iy37+7mVojJXyB/68zs6bAIqBzYWEhPXr0YO3ateQkVsLOycmhrKysxvcWFxdTXFwMwLp164IuVcS7+24/mVS7dnXP8GcGRUVw6qnh1SayC/MXySEcyKxVfn5+xZgxY+jVqxcbN27c+Vrr1q3rbQfPy8tj4cKFAVcpWe+DD+CEE+Dss+Evf9EUrZIOav1PGNpshM65jfn5+cycOZMOHTpQWloKQGlpKe3btw+rDJHaOQdXXOHn177vPoW3pL2ge6G0M7NWie9bvvTSSxx//PH079+fkpISAEpKShgwYECQZYgk57HHYM4c3yxy2GFRVyNSr6DbwHOAkkQ7eJPCwkL69etHz549GTRoEOPHj6djx45Mnjw54DJEEt59F371K/jiiz1fe/ttOOMMfxUuEgOhtYEn7PXB1AYujVZdDT17wvvvQ17enq/vv79f1b1z5/BrE6ldrW156vsk2WPsWL96+8SJ8NOfRl2NSKNpSTXJDh9/7LsG9u0LF1wQdTUiKaEAl8znHFx1lf/+gQfUu0QyhppQJHOMGeNn/dt9zUnn/E3Le+6BI4+MpDSRICjAJTOsWAFDh8Jpp/mVb3Z3xBF+3hKRDKIAl/jbvh0uv9z3Ipk6FTQwTLKEAlzir7gYXn8dJkxQeEtWUYBLPFRW+pVvdldRATfeCP/5n/Czn4Vfl0iEFOCS/ioroVcvmD+/5tdbtoQHH1TvEsk6CnBJf3fd5cN7xIiae5Gccgocc0zoZYlETQEu6W3lSh/cAwdCYkk+EfE0kEfS147pXffd1/fxFpFv0BW4RGPpUj+8vS6LFsHLL8O4cTWv9i6S5RTgEr65c+Gss5LbNz8fLrssyGpEYksBLuHauhWGDIGjj4ann4Ym9bTinXhi/fuIZCkFuITrllv8jcnZs2uek1tEkqZLGwnP22/D7bfD4MF+4I2INIquwCU4X30FjzwCGzb455MnQ+vWvl+3iDSaAlyCc8st8Mc/fv28RQu/Gk7bttHVJJJBFOASjGXL/OruF14Ijz7qt5lBM/2XE0kV/TRJ6lVX+65/Bx4Io0ZB8+ZRVySSkZK+iWneRWb2u8TzjmZWw8z5kvUeeADefNOHd7t2UVcjkrEacgV+P7Ad6A38AdgMPAOcFkBdku6WL4c//cnfqNzdCy9AYSFcdFH4dYlkkYYEeA/n3ClmtgTAOVdhZvsEVJeks6++gp/8BD76CA4/fM/XTzrJD3/X9K4igWpIgFeaWVPAAZhZO/wVuWSbO+6Ad96B6dPh3HOjrkYkazVkIM9o4C9AezO7Ffg7cFsgVUn6eu893zVw0CCFt0jEkr4Cd85NNLNFQAFgwH8551YEVpmkn+3b/TwmLVvCvfdGXY1I1ks6wM1sPDDGOXffLttGOOdGBFGYRGj0aN+TZHdffQUffggPPwyHHBJ+XSLyDeacS25Hs9XAeuBu59xjiW2LnXOnNOB4yR2sBnl5eSxcuHBv3y7JWrwYTjvNL1N29NF7vt61K9x0k25QioSn1h+2htzELAPygYlm1gO4rq4PlhiqqvIDcNq397MFtmoVdUUiUoeG3MQ059wm59y5wDrgFeCgYMqSSIwaBUuWwNixCm+RGGjIFfj0Hd8450aY2ULg+tSXJIFyzl9p727VKr9o8IABfgFhEUl7SbeBp4jawKP05Zdw5pkwf37Nrx9wAKxYAYcdFm5dIlKXvW8DN7O/O+d6mdlmvhnABjjn3IEpKFDCcMstPryHDq25ieT731d4i8RIvQHunOuVeDwg+HIkMG+/DSNH+tVw7rgj6mpEJAUaMhvhMWbWIvF9vplda2atAqtMUqe6Gi6/XKvhiGSYhtzEfAbIM7POwHj8Tc0ngR8EUZg00r/+BZ995r9/7jnfdKLVcEQySkMCfLtzrsrMfgjc45wbs2NmQkkzy5b5Fd+//PLrbX37wgUXRFeTiKRcQ2cjvAAYDOyYxUhLraSb7dt9c8l++/kr7mbNoGlT6N1boydFMkxDBvJcDPQEbnXOrTKzo4An6nqDmR1hZi+b2QozW3ZvYgKk8vJyCgsLyc3NpbCwkIqKir0+AdnNAw/AvHl+UM7AgdC/P5xzjp+ASkQySsr6gZvZM865H+22LQfIcc4tNrMDcnNzNz377LNMmDCBNm3aMGzYMIqKiqioqGDkyJF1fr76gSfh00/h29+Gnj1h1ixdcYtkhpTMhVKfPWY+cs6VAqWJ7zcPGDCANWvWMG3aNObOnQvA4MGDyc/PrzfApQbOwYwZUFbmnz/xhO9x8uCDCm+RLJDKAK/zUt7MOh1xxBH06NGDtWvXkpOTA0BOTg5lOwJoN8XFxRQXFwOwbt26FJaaIcaNg6uu+ua20aNrnkVQRDJOKptQap1a1sz2B1555plnThk4cCCtWrVi48aNO19v3bp1ve3gakLZzZo1vrkkLw/Gj/fbWrSAxC9GEckYtf453ZCbmHt1EDNrju9DPnFgYpKkDh06UFpaCkBpaSnt27dPYRlZ4pe/9AssjBsHnTr5L4W3SFZJZYDfuPsGMzP8oJ8Vzrm7d2zv378/JSUlAJSUlDBgwIAUlpEFpk6Fv/wFRoyAzp2jrkZEIlJvE4qZ/YM62redc93qeG8v4DXgH8D27t27d7/tttvo0aMHgwYN4pNPPqFjx45MnjyZNm3a1FlH1jWhlJbCmDHfHIyzw6RJfkmz+fOhubrii2S4RvVC6Zd4vCbx+Hji8ULg33W90Tn3990OvvMXwZw5c5I4dJZyDi68EF55xQ/I2V2bNr7dW+EtktWSmY3wYwAzO8M5d8YuLw0zs9eBPwRVXNZ69FF4+WXfHfCKK6KuRkTSVEPawPdLNIkAYGanAzVcHkqjrF3r5+v+3vf8kHgRkVo0pB/4pcAjZrZjHcyNwCUpryjbXXcdbNkCxcXQJJX3mEUk0yQd4M65RUB3MzsQf/Pz8+DKyiIvvgh33uknoaqshFdfhT/8AY4/PurKRCTNJbOk2kXOuSfM7PrdtgOwa/dAaaCyMj/Fa8uWvh83wMUXw4179MgUEdlDMlfgO9q5taRaqv3617B5M7z2mh9VKSLSAMn0QhlnZk2BTc65USHUlB1mzIAnn4Sbb1Z4i8heSeoumXOuGugfcC3Z44sv/CRUXbrA8OFRVyMiMdWQXihvmNlY4Glgy46NzrnFKa8qEw0b5vt3A2zbBp9/7ptOWrSIti4Ria2GBPjpicffJx4NP7Kyd0orykQvvggjR0JhIRxzjN/Wq5f/EhHZSw0J8Ofwgb1jaLwDNpnZSc65t1JdWMbYsgWuvBKOPRamT4d99426IhHJEA0J8FOBPGA6PsTPARYAV5jZZOfc7QHUF38jRsCqVX5eE4W3iKRQQwK8LXCKc+4LADO7GZgC/AewCFCA727xYrj7bj8k/j/+I+pqRCTDNGSsdkfgq12eVwJHOue2AttSWlUmqKqCyy6D9u3hdv1uE5HUa8gV+JPAm2Y2LfH8XGCSme0HLE95ZXE3ahQsWQKTJ0OrVlFXIyIZqEFrYprZqUAvfBv4351zDV1hYa8X4IzVgg4ffggnnADf/75fOUcrxIvI3mvUgg47JSa0WtTocjKZc34O72bN4L77FN4iEpgGBbjU4pNP4L33/PcLFsBLL/nwPuywaOsSkYymAG+sjz7yzSVbtny9rVcv3/dbRCRACvDGcM4HtRnMmvX1+pV5eVqMQUQCpwBvjIkTfXCPHu1vWIqIhEiXiXtr/Xo/n3ePHnD11VFXIyJZSFfgyXIOpk71bd4As2fDxo3w0EPQtGmUlYlIllKAJ+vpp/3yZ7v605/gxBOjqUdEsp4CPBkbNsC118Jpp/kr7yZN/NeOm5YiIhFQgCdj6FAoL/fhfdBBUVcjIgLoJmb95syBCRPgN7+B7t2jrkZEZKcGzYWSAvGaC2XrVt/GbQZvvw0tW4Z7fBGRVM2FknV+/3v44AP4298U3iKSdtSEUpslS+DOO+GSS+Css6KuRkRkDwrwmlRV+VV02raFO+6IuhoRkRqpCaUmo0fDokXw1FPQpk3U1YiI1EgBDrBiBVx1FWza5J8vXw79+sGgQdHWJSJSBwV4dTX8/OewcqWfBhaga1cYOVKLMYhIWlOA33cfzJ8PTz6551B5EZE0lt03MT/+GP7nf6BvXzj//KirERFpkOwNcOe+ngb2gQfUXCIisZO9TShPPw0vvAD33ANHHhl1NSIiDZadV+C7zi74i19EXY2IyF7JzivwoUOhosKvHq/FGEQkpgK9AjezR8yszMze2bGtvLycwsJCcnNzKSwspKKiIsgS9rRjdsEbboBu3cI9tohICgXdhDIB6LPrhqKiIgoKCli5ciUFBQUUFRUFXMIu/v1vuOIKyM2Fm24K77giIgEINMCdc68C5btumzZtGoMHDwZg8ODBPPvss0GW8E07ZhcsLoZ99w3vuCIiAQi9DXzt2rXk5OQAkJOTQ1lZWa37FhcXU1xcDMC6desad+AlS+Cuu+DSSyE/v3GfJSKSBgJf0MHMOgHPOedOAFyrVq3YuHHjztdbt26dVDt4oxZ0qKqC734XVq/28560br13nyMiEr5aB6mE3o2wQ4cOlJaWAlBaWkr79u2DP+iO2QXHjFF4i0jGCD3A+/fvT0lJCQAlJSUMGDAg2AOuWuVvWJ57Lpx3XrDHEhEJUdDdCCcB84DjzGz1+PHjGTZsGLNnzyY3N5fZs2czbNiw4Apwzk8T26SJn7RKw+VFJIMEehPTObf79H4OYM6cOUEe9msTJ8KsWb7p5IgjwjmmiEhIMnco/fr18OtfQ8+e/ipcRCTDZG6AX389fP45PPSQhsuLSEbKzAB/8UV4/HEYNsyvriMikoEyL8C3bPHD5Y87zi/WICKSoTJvNsKbb4aPPoJXX9VweRHJaJl1Bb5oEYwa5a/Av/e9qKsREQlU5gR4ZSVcdhl06OBXlBcRyXCZ04QyahS89RZMnQoHHRR1NSIigcuMK/APPvBt3z/8of8SEckC8Q9w53yb9z77wNixUVcjIhKa+DehlJT4ZdIeeAAOPTTqakREQhPvK/C1a/2Iy169YMiQqKsREQlVvAN8+HA/cKe42M84KCKSReKbelu3wtNPw+DB0KVL1NWIiIQuvgE+e7ZfZV6LNIhIlopvgE+dCq1aaYFiEcla8Qzwykr461+hXz/ffVBEJAvFM8BffRXKy2HgwKgrERGJTDwDfOpUaNkSzj476kpERCITvwDfvh2efRb69IFvfSvqakREIhO/AJ8/Hz77TM0nIpL14hfgU6dCs2ZwzjlRVyIiEqn4BfjMmXDmmdC6ddSViIhEKl4BXlUF774Lp54adSUiIpGLV4B/+KHvA66h8yIiMQvwFSv8owJcRCRmAf7uu/7x+OOjrUNEJA3EK8BXrICcHK15KSJCHANczSciIoACXEQktuIT4JWVsHmzAlxEJCE+Ab51q39UgIuIAHEK8C+/9I8KcBERIG4BfuCBcMghUVciIpIW4hXgXbqAWdSViIikhfgE+Nataj4REdlFPAJ840Y/kZUCXERkp3gEuOZAERHZQ7wCXHOgiIjsFFmAz5w5k+OOO47OnTtTVFRU984rVvibl0cdFU5xIiIxEEmAV1dXc8011zBjxgyWL1/OpEmTWL58ee1vePddaNHCL6UmIiJARAE+f/58OnfuzNFHH80+++zD+eefz7Rp02p/w4oV0LJleAWKiMSAOedCO1ifPn3c+vXrqaioYNOmTRx55JEAbNiwgS1bttCxY8dv7L9u3TrWr1/PwZWVlFZX0+3kk0OrNR2sW7eOdu3aRV1GqLLtnLPtfEHn3FCLFi2a5ZzrU9NroQY44AAmT57MrFmzePjhhwF4/PHHmT9/PmPGjKn1jfvttx9btmwJp8o0kZeXx8KFC6MuI1TZds7Zdr6gc94LtY5ejKQJ5fDDD+fTTz/d+Xz16tUceuihUZQiIhJbkQT4aaedxsqVK1m1ahVfffUVTz31FP3794+iFBGR2IqkW0ezZs0YO3YsZ599NtXV1VxyySV07dq1zvccfPDBIVWXPoYMGRJ1CaHLtnPOtvMFnXMqRdIGvjeysd1MRIR0awMXEZHGU4CLiMRUWgV4fcPrnXNce+21dO7cmW7durF48eIIqkyt+s554sSJdOvWjW7dunH66aezdOnSCKpMrWSnUViwYAFNmzZlypQpIVYXjGTOee7cuZx00kl07dqVM888M+QKU6++c/78888599xz6d69O127duXRRx+NoMrUueSSS2jfvj0nnHBCja8Hkl/OuTC/alVVVeWOPvpo98EHH7ht27a5bt26uWXLlu18/dRTT3XPP/+869Onj9u+fbubN2+e+853vlPXR6a9+s7ZOedef/11V15e7pxz7oUXXsiKc96x31lnneX69u3rJk+eHEGlqZPMOVdUVLguXbq4jz/+2Dnn3Nq1a6MoNWWSOedbb73V3XDDDc4558rKylzr1q3dtm3boig3JV555RW3aNEi17Vr1xpfb0R+1ZqpYQd47YVAT2DWLs+HA8N3eT4TGAdcsMu294CcqGsP6pxr2L81sCbqusM4Z+BXwDXABOC8qOsO+pyBq4Fboq415HMeDtyPv0l3FPBPoEnUtTfyvDsB79TyWsrzK52aUA4DPt3l+erENgCcH0pa5z4x1NDzuRSYEWhFwav3nM3sMOCHwIMh1hWkZP6djwVam9lcM1tkZj8LrbpgJHPOY4EuwGfAP4DrnHPbwykvEinPr3Sa3q+mrjK7dztMZp84Sfp8zOwsfID3CrSi4CVzzvcANzrnqi0z1kBN5pybAacCBUBLYJ6Zvemcez/o4gKSzDmfDbwF9AaOAWab2WvOuU0B1xaVlOdXOgX4auCIXZ4fjv/N3NB94iSp8zGzbsDDQF/n3IaQagtKMuecBzyVCO+DgR+YWZVz7tlQKky9ZP9vr3fObQG2mNmrQHcgrgGezDlfDBQ5357wTzNbBRwPzA+nxNClPL/SqQllAZBrZkeZ2T7A+cD03faZDvzMvO8CnzvnSsMuNIXqPWcz6whMBf5fjK/GdlXvOTvnjnLOdXLOdQKmAFfHOLwhuf/b04DvmVkzM/sW0ANYEXKdqZTMOX+C/4sDM+sAHAd8GGqV4Up5fqXNFbhzrsrMfgHMApoCjzjnlpnZlYnXHwReAH6Av9nxb/xv8NhK8px/B7QF7k9ckVY55/KiqrmxkjznjJLMOTvnVpjZTOBtYDvwsHPuneiqbpwk/53/CEwws3/gmxdudM6tj6zoRjKzSUA+cLCZrQZuBppDcPkV9lB6ERFJkXRqQhERkQZQgIuIxJQCXEQkphTgIiIxpQAXEQmImT1iZmVmVm+PIjM70szmmNnbiRG5h9f3HgW4ZB0z62RmP23kZ/wq0V9bpC4TgBpXlK/BncBjzrluwB+AP9X3BgW4ZKNOQKMCHD/ZlgJc6uScexUo33WbmR1jZjMTc968ZmbHJ176NjAn8f3LwID6Pl8BLhnDzC4ys/lm9paZjTOzHok/R/c1s/3MbJmZnQAU4Uc9vmVmvzazpmZ2h5ktSOx/ReLz8hN/yk4xs3fNbGJiFN21wKHAy2b2cpTnLLFUDPzSOXcqMBQ/IyPAUuBHie9/CBxgZm3r+qC0GYkp0hhm1gX4CXCGc67SzO7HD82eDtyCnyDqCefcO2Y2DBjqnOuXeO8Q/LDm08ysBfC6mb2Y+OiTga74OSteT3z+aDO7HjgrziMHJXxmtj9wOjB5l4naWiQehwJjzeznwKvAGqCqrs9TgEumKMDP5rcg8YPREijDtyUuAL4Erq3lvd8HupnZeYnnBwG5wFfAfOfcagAzewvf/PL3QM5AskETYKNz7qTdX3DOfQYMhJ1B/yPn3Od1fZgCXDKFASXOueHf2Gh2CLA/fk6KfYEttbz3l865Wbu9Nx/YtsumavQzI43gnNtkZqvM7MfOucnmrza6OeeWmtnBQHliTvThwCP1fZ7awCVTzAHOM7P2AGbWxsyOxLc33gRMBEYm9t0MHLDLe2cBV5lZ88R7jzWz/eo53u6fIbKHxARX84DjzGy1mV0KXAhcamZLgWV8fbMyH3jPzN4HOgC31vf5upqQjOCcW25mvwVeNLMmQCV+itYq59yTZtYUeMPMegOvAVWJH6AJwL34ppHFiSuidcB/1XPIYmCGmZU6584K4pwk/pxzF9Ty0h5dC51zU/DTJydNsxGKiMSUmlBERGJKAS4iElMKcBGRmFKAi4jElAJcRCSmFOAiIjGlABcRian/DzYf8KmpShy7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "Extent_area = np.linspace(1000,1000000000,100)\n",
    "grid_size = calcMingrid(Extent_area)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.xlabel(\"extent\")\n",
    "plt.ylabel(\"grid_size\")\n",
    "\n",
    "# plot the function\n",
    "plt.plot(Extent_area, grid_size, 'r')\n",
    "\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "b3ab0ba2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calcMingrid(153531.69540990253)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64e916c2",
   "metadata": {},
   "source": [
    "We can say that for square grids, we can use our function to aprox the ratio between extent and min grid size because smaller grids definitions are more accurate and bigger grid size in relation with the extent. Also our function is eager that the real calculation done by posgis so we can set a smaller limit in our function. Also as a minimum guide to choose grid size we will take into account a relation between km/pixel and the extent it represents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84814a63",
   "metadata": {},
   "source": [
    "# Queries to generate grids on the fly:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbbfd96a",
   "metadata": {},
   "source": [
    "`ST_SquareGrid` or `ST_HexagonGrid`\n",
    "```sql\n",
    "select geom as the_geom_webmercator from (SELECT (ST_SquareGrid(4000, ST_Transform(a.the_geom, 3857))).* \n",
    "FROM admin_regions a  \n",
    "WHERE gid_0 = 'ESP' and gid_1 is null and gid_2 is null) grid;\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96a4d2d1",
   "metadata": {},
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b4932b4",
   "metadata": {},
   "source": [
    "```sql\n",
    "INSERT INTO planning_units_geom   (the_geom, type, size)    \n",
    "select st_transform(geom, 4326) as the_geom,'square' as type, 25 as size from \n",
    "(SELECT (ST_SquareGrid(25000, ST_Transform(ST_GeomFromGeoJSON('{\"type\":\"Polygon\",\n",
    "        \"coordinates\":[[[-78.046875,-36.3151251474805],[118.47656249999999,-36.3151251474805],[118.47656249999999,75.23066741281573],[-78.046875,75.23066741281573],[-78.046875,-36.3151251474805]]]}'), 3857))).* \n",
    " ) grid\n",
    " ON CONFLICT ON CONSTRAINT planning_units_geom_the_geom_type_key DO NOTHING;\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
